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SUMMARY 


This  research  work  developed  a  systematic  methodology  for  the  design  of  op¬ 
timal  process  parameters  for  forging  processes.  The  methodology  uses  a  multi¬ 
disciplinary  approach  involving  advanced  modeling  methods,  finite  element  tech¬ 
niques,  nonlinear  mechanics  and  modern  control  theory,  for  maintaining  the  forg¬ 
ing  system  variables  within  ‘favorable’  processing  windows.  Appendix  A  gives 
an  overview  of  the  project  goals  and  tasks,  and  also  describes  in  brief  the  tasks 
accomplished  in  this  work  and  their  organization  in  this  report.  Chapter  1  in¬ 
troduces  the  objectives  and  goals  of  this  research  work,  and  the  design/ analysis 
methods  used  to  achieve  these  goals.  It  also  includes  a  section  surveying  the  ear¬ 
lier  work  done  in  this  area.  Chapter  2  gives  a  brief  description  of  the  nonlinear 
rigid  viscoplastic  finite  element  method,  and  presents  the  mathematical  deriva¬ 
tion  of  the  deformation  state  space  model  from  the  finite  element  equations.  In 
addition,  this  chapter  also  describes  the  finite  element  based  condensation  pro¬ 
cedure  used  to  reduce  the  number  of  degrees  of  freedom  of  the  deformation  state 
space  model.  The  procedure  for  development  of  the  thermal  state  space  mod¬ 
el  from  the  finite  element  governing  equations  is  described  in  Chapter  3  of  the 
report.  This  chapter  also  presents  the  procedure  for  condensation  of  the  ther¬ 
mal  state  space  model  to  reduce  the  order  of  the  state  space  system.  Chapter 
4  describes  the  development  of  the  coupled  thermomechamcal  state  space  mod¬ 
el,  and  the  construction  of  the  control  output  matrix.  Chapter  5  explains  the 
design  strategy  used  for  optimizing  the  process  parameters  based  on  the  linear 
quadratic  regulator  (LQR)  theory.  This  chapter  gives  a  detailed  description  of 
the  procedure  used  in  solving  the  state  space  equations,  and  also  presents  the 
methodology  used  in  designing  the  initial  die  temperature.  The  logic  used  in 
selecting  weighting  matrices  for  the  LQR  control  scheme  is  also  described  in 
Chapter  5.  The  methodology  developed  was  tested  and  validated  using  a  variety 
of  forging  simulations  under  different  operating  conditions.  These  test  cases  and 
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numerical  examples  are  presented  in  Chapter  6  of  this  report.  While  dealing 
with  large  scale  finite  element  models  involving  a  large  number  of  degrees  of  free¬ 
dom,  generally,  numerical  difficulties  are  encountered  during  the  implementation 

•  process.  To  avoid  such  situations  and  to  reduce  the  time  and  effort  during  com¬ 
putation,  reduced  order  models  need  to  be  used.  Chapter  7  describes  in  detail 
the  model  reduction  methods  studied  and  utilized  during  the  course  of  this  work. 

*  This  chapter  also  presents  the  comparative  evaluation  of  several  order  reduction 
methods  and  the  basis  for  selection  of  the  Balanced  Model  Reduction  (BMR)  as 
suitable  for  use  in  this  work.  The  control  law  and  weighting  matrix  selection  had 
to  be  modified  for  use  with  the  reduced  order  models.  The  procedure  for  doing 
this,  and  a  few  case  studies  to  validate  the  conclusions  drawn  are  also  presented 
in  Chapter  7.  Chapter  8  presents  more  complicated  and  large-scale  simulation 
examples  to  demonstrate  the  effectiveness  of  the  model  reduction  methodology 
developed.  After  the  control-based  methodology  was  developed,  approximations 
were  made  and  introduced  into  the  program  to  reduce  the  computational  time 
involved  in  the  process.  Chapter  9  presents  the  results  of  these  studies  and  de¬ 
scribes  the  effect  of  using  one  state  space  model  for  the  entire  simulation,  and 
the  effect  of  using  smoothened  (approximate)  die  velocity  schedules  during  the 
simulation. 

To  summarize  the  entire  procedure,  the  nonlinear  rigid  viscoplastic  finite 
element  method  is  used  for  deformation  and  thermal  analyses.  A  coupled  state 
space  model  is  then  built  to  represent  the  deformation  and  thermal  behavior  of 
the  material,  with  nodal  velocities  and  nodal  temperatures  as  the  state  variables. 
A  finite  element  based  condensation  technique  is  used  for  reducing  the  order 
of  the  system  to  facilitate  numerical  analysis.  Sophisticated  model  reduction 
techniques  are  used  to  further  reduce  the  order  of  the  state  space  system.  The 
desired  (favorable)  forging  conditions  are  obtained  by  imposing  constraints  on 

%  effective  strain-rate  and  temperature  variation  within  the  deforming  material. 

The  linear  quadratic  regulator  (LQR)  theory  with  finite  time  control  is  used  in 
designing  the  ram  velocity  and  initial  die  temperature. 
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CHAPTER  1 


Introduction 

« 


1.1  Need  for  process  control  in  metal  forming 

In  today’s  metal  forming  industry  new  materials  are  being  widely  introduced 
and  used  for  an  ever-broadening  scope  of  applications.  The  increasing  demand  for 
new,  high  performance  materials  has  augmented  the  need  for  superior  processing 
methods.  Computer-aided  design  techniques  have  generally  been  found  to  be  the 
most  effective  and  efficient  way  to  meet  this  challenge  and  are  used  widely.  In 
this  work,  the  nonlinear  rigid  viscoplastic  finite  element  method  has  been  used 
for  simulation  and  analysis  of  forging  processes.  Even  though  this  work  focuses 
on  forging,  the  methodology  featured  in  this  report  is  general  purpose  and  is 
applicable  to  any  of  the  unit  forming  processes. 

During  manufacturing  the  mechanical  and  service  properties  of  the  final  prod¬ 
uct  are  dependent  to  a  large  extent  on  the  prior  processing  history.  As  such,  one 
of  the  most  important  tasks  is  the  selection  and  control  of  critical  process  pa¬ 
rameters  that  would  ensure  required  part  quality  along  with  specific  mechanical 
'  and  physical  characteristics.  Therefore,  there  is  a  need  to  develop  an  optimal 

processing  strategy  that  would  result  in  defect-free,  high  quality  products  on  a 
»  repeatable  basis. 
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Forging  is  a  thermomechanical  plastic  deformation  process  wherein  the  work 
piece  is  deformed  from  a  relatively  simple  geometry  to  a  predetermined  complex 
shape  by  the  application  of  compressive  forces.  Forging  differs  from  other  shap¬ 
ing  methods  in  that  the  flow  of  metal  is  intended  to  produce  specific  material 
properties  in  the  final  product.  The  forging  operation  can  be  visualized  as  a 
system  with  a  large  number  of  interacting  factors  like  starting  billet  shape,  inter¬ 
face  frictional  conditions,  temperature  of  workpiece  and  dies,  velocity  of  the  die, 
and  geometry  of  the  final  part.  These  parameters  strongly  influence  the  thermal 
and  flow  behavior  of  the  deforming  material  and  have  a  direct  impact  on  the 
spatial  and  temporal  distribution  of  field  variables  like  strain-rate,  total  effective 
strain,  and  nodal  temperatures.  The  properties  and  integrity  of  the  final  formed 
product  are,  in  turn,  functions  of  the  field  components  of  the  metal  forming 
system.  It  is  thus  necessary  to  monitor  and  control  the  essential  field  variables 
like  strain,  strain-rate,  and  temperature  to  obtain  the  desired  microstructure  and 
service  properties  in  the  final  product.  Also,  non-conventional  and  difficult-to- 
process  materials  generally  have  a  narrow  range  of  processing  conditions  in  the 
strain,  strain-rate,  and  temperature  space  where  they  can  be  processed  success¬ 
fully  without  any  formation  of  defects.  To  maintain  the  field  variables  within 
these  ‘favorable’  processing  regions,  optimal  design  of  process  parameters  such 
as  die  velocity  and  initial  die/billet  temperature  must  be  carried  out. 
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1.2  Process  parameters  and  physical/microstructural  properties 


Generally,  metal  forming  processes  may  be  divided  into  two  categories:  mas¬ 
sive  forming  and  sheet  metal  forming.  The  massive  forming  processes  include 
operations  such  as  forging,  extrusion,  and  rolling.  Although  this  work  is  generic 
in  nature  and  is  applicable  to  all  the  unit  forming  processes,  the  main  focus  is 
on  isothermal / nonisothermal  forging.  In  these  processes  there  is  a  strong  rela¬ 
tionship  between  macroscopic  (physical)  properties,  microstructure,  and  process 
parameters.  This  is  strongly  supported  by  earlier  work  done  in  this  area.  For  in¬ 
stance,  in  the  design  of  aircraft  engine  disks,  attempts  have  been  made  to  obtain 
desired  properties  in  the  final  product  by  generating  different  microstructures 
in  different  regions  of  the  disk  [1].  Though  detailed  relationships  between  mi¬ 
crostructure  and  physical  properties  have  not  yet  been  fully  established,  it  has 
been  shown  through  experiments  that  physical  properties  depend  to  a  large  ex¬ 
tent  on  the  microstructure  of  the  material,  and  if  the  workpiece  is  not  processed  in 
the  right  manner  and  under  appropriate  operating  conditions,  defective  products 
and/or  unwanted  physical  characteristics  in  the  final  product  may  result. 

In  metal  forming  processes,  it  has  generally  been  found  that  process  param¬ 
eters  play  an  important  role  in  the  forming  of  a  specific  microstructure  in  the 
workpiece.  For  example,  Devadas  et.  al.  [2],  modeled  the  microstructure  and 
mechanical  properties  of  steel  during  hot  rolling.  These  models  have  been  used 
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to  predict  grain  size  and  characterize  static  and  dynamic  recrystallization  dur¬ 
ing  the  rolling  process.  Most  of  these  models  have  been  found  to  be  sensitive  to 
change  in  process  parameters  such  as  strain,  strain-rate,  and  temperature.  Sastry 
et.  al.  [3]  studied  the  relationship  between  metallurgical  properties  and  process 
parameters  during  the  superplastic  forming  of  titanium  alloys.  This  work  showed 
that  the  percentage  of  equilibrium  a  phase  in  the  microstructure  is  related  to 
strain-rate,  strain-rate  sensitivity,  and  temperature.  Dadras  and  Thomas  [4]  con¬ 
ducted  experiments  to  study  the  deformation  behavior  of  Ti-6242  (Ti-6Al-2Sn- 
4Zr-2Mo-0.lSi)  alloy  during  the  upset  forging  of  specimens  starting  with  (ol  +  P) 
or  p  microstructure.  They  found  that  the  volume  percentage  of  primary  a  mi- 
crostructure  in  the  specimen  is  influenced  by  the  deformation  temperature.  This 
further  illustrated  and  emphasized  the  importance  of  temperature  effects  during 
deformation.  Recently,  Cohen  and  Durham  [5]  correlated  the  effect  of  change  in 
strain-rate  to  the  final  microstructure  during  hot  working  processes.  The  depen¬ 
dence  of  the  resulting  microstructure  on  various  temperatures  and  strain-rates 
was  analyzed  while  carrying  out  compression  tests  on  a  carbon  steel  material. 
Seetharaman  et.  al.  [6]  analyzed  the  effect  of  strain,  strain-rate,  and  temper¬ 
ature  on  the  microstructure  of  a  gamma  Ti-Al  alloy  during  extrusion.  During 
this  work,  the  relationship  between  microstructural  parameters  (grain  size  and 
grain  distribution)  and  process  parameters  (effective  strain-rate,  effective  strain, 
and  flow  stress)  was  obtained  using  experimentally  collected  data.  This  further 
proved  that  both  strain-rate  and  temperature  have  a  strong  influence  on  the 
microstructure  and  properties  of  a  given  material. 
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Besides  microstructure,  workability  is  also  an  important  characteristic  in  an¬ 
alyzing  metal  forming  processes.  Workability  is  the  capacity  of  a  material  to 
deform  without  failure.  It  depends  on:  (1)  process  parameters  (such  as  temper¬ 
ature,  strain-rates,  stresses,  and  strain  history),  and  (2)  material  variables  (such 
as  composition  and  initial  microstructure).  The  current  work  proposes  that  it  is 
crucial  to  select  proper  processing  conditions  and  process  parameters  to  achieve 
the  desired  microstructure  and  workability  level  in  the  deforming  workpiece.  For 
a  given  billet  material  and  geometry  the  following  major  process  parameters  have 
to  be  determined  to  obtain  optimal  forging  conditions. 


•  Die  material 


•  Initial  die/ workpiece  temperatures 


•  Ram  speed 


•  Lubricant 


•  Preform  geometry  and  position 


•  Die  geometry 


•  Type  and  size  of  the  forging  press 
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In  the  past,  the  design  of  process  parameters  for  metal  forming  operations 
was  done  on  a  trial  and  error  basis,  and  was  dependent  to  a  large  extent  on  the 
experience  of  the  die  designers,  metallurgists  and  process  engineers.  This  method 
was  expensive,  time  consuming,  and  not  very  reliable.  With  the  development  of 
sophisticated  numerical  techniques,  computer  aided  design  and  analysis  methods 
have  become  very  popular  and  are  in  use  widely. 

1.3  Metal  forming  simulation  using  the  finite  element  method  (FEM) 

With  advancement  in  computer  technology  and  development  of  sophisticated 
processing  methods,  the  use  of  computer-aided  techniques  for  process  simulation 
and  process  design  has  increased  considerably.  Some  of  the  well  known  metal 
forming  simulation  and  analysis  methods  are  the  slab  method,  the  upper-bound 
method,  the  slip-line  field  method,  and  the  finite  element  analysis  (FEA)  method. 

The  slab  method  is  based  on  the  equilibrium  of  a  slab  of  the  deforming  body 
and  assumes  a  simplified  stress  distribution  along  the  slab.  Biswas  and  Rooks  [8] 
evaluated  the  loads  and  stresses  for  metal  forming  processes  based  on  this  method 
using  an  approach  in  which  the  various  deformation  stages  are  uncoupled  and 
analyzed  separately.  Lui  and  Das  [9]  also  used  the  slab  method  for  evaluating 
the  loads  and  stresses  in  axisymmetric  forgings.  Although  this  method  is  quick 
and  gives  reasonable  results  [10,11],  its  main  shortcoming  is  that  it  is  restricted 
to  the  evaluation  of  loads  and  stresses  for  simple  geometries  only. 


* 
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The  upper  bound  method  [12,13]  assumes  a  velocity  field  describing  metal 
flow  during  the  forming  operation.  Based  on  this  velocity  field,  the  total  forming 
energy  and  the  forming  load  are  then  calculated.  By  including  one  or  more  pa¬ 
rameters  in  the  considered  velocity  field,  it  is  possible  to  determine  and  optimize 
the  upper-bound  velocity  field.  The  optimal  process  parameters  are  determined 
by  minimizing  the  total  forming  energy  with  respect  to  these  parameters  [14]. 
This  approach  has  earlier  been  used  in  the  analysis  of  axisymmetric  forgings, 
plate  rolling,  and  extrusion  processes  [15-19]. 

Due  to  the  rapid  development  of  computers  and  numerical  algorithms,  the 
finite  element  method  (FEM)  has  become  very  popular  for  simulating  and  ana¬ 
lyzing  metal  forming  problems.  Most  of  the  FE  analyses  in  this  field  are  based 
on  the  rigid-plastic  and  rigid-viscoplastic  theories.  For  rigid-plastic  materials,  it 
is  assumed  that  flow  stress  is  a  function  of  strain,  strain-rate,  and  temperature, 
and  that  the  elastic  response  of  the  material  is  negligible.  The  rigid- viscoplastic 
theory  was  generalized  by  Zienkiewicz  et.  al.  [20,21]  and  is  capable  of  model¬ 
ing  hot,  rate-dependent  processes.  This  generalization  provides  the  theory  for 
analyzing  the  deformation  of  Ti-alloys  which  are  strain-rate  sensitive  materials. 

Lee  and  Kobayashi  [22,23]  developed  the  rigid-plastic  finite  element  method 
using  variational  principles  for  a  material  obeying  Von  Mises’  yield  criterion  with 
isotropic  kinematic  hardening.  In  the  early  80’s,  Oh  et.  al.  [24,25]  refined  the 
rigid- viscoplastic  formulation  to  solve  a  wide  variety  of  problems  using  the  fi¬ 
nite  element  method.  These  efforts  resulted  in  the  development  of  a  generalized 
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computer  program  ALPID  [26]  (Analysis  of  Large  Plastic  Incremental  Deforma¬ 
tion),  which  has  the  capability  to  perform  a  wide  range  of  2-D  metal  forming 
simulations. 

A  variety  of  metal  forming  problems  have  been  solved  using  the  finite  ele¬ 
ment  analysis  method.  Reference  27  gives  a  broad  and  detailed  coverage  of  the 
application  of  FEM  in  metal  forming  research.  It  contains  plane  strain  problems 
(such  as  bulk  forging,  sheet  rolling,  plate  bending  and  side  pressing),  axisymmet- 
ric  forgings  (such  as  disk  forging  and  ring  compression),  steady-state  processes 
such  as  extrusion  and  drawing;  sheet  metal  forming  operations,  and  the  forging 
of  porous  metals.  Duggirala  [28]  also  used  ALPID  in  the  analysis  of  flashless  ring 
gear  blanks  and  axle  shafts. 

Besides  the  rigid-plastic  and  rigid- viscoplastic  finite  element  analysis,  some 
researchers  also  used  other  plastic  theories.  Dexter  [29]  studied  the  mechanisms 
involved  in  forging  using  the  elastic- viscoplastic  finite  element  analysis.  The 
stress  and  strain  values,  especially  at  the  die- workpiece  interface,  and  the  critical 
loads  for  the  potential  failure  of  forging  die  can  be  obtained  by  this  method  of 
analysis.  The  updated  Lagrangian  Jaumann  Formulation  of  FEM,  which  includes 
elastic  deformation  and  tends  to  be  elastic-plastic  or  elastic- viscoplastic  in  nature, 
has  also  been  used  in  metal  forming  [30].  It  has  successfully  been  applied  to  the 
problems  of  extrusion,  drawing,  rolling  and  sheet  metal  forming. 

Finite  element  analysis  and  simulation  has  been  further  enhanced  and  devel¬ 
oped  by  several  researchers  to  handle  nonisothermal  processing  conditions.  Wu 
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and  Oh  developed  the  program  ALPIDT  [31]  which  interfaces  ALPID  with  a 
thermal  analysis  module,  and  can  perform  coupled  thermo- viscoplastic  deforma¬ 
tion  analyses.  Oh  et.  al.  [32]  used  ALPIDT  to  analyze  the  hot  die  forging  of 
Ti6242  alloys.  They  found  that  the  distortion  in  the  FE  mesh  was  more  severe 
during  hot-die  forging  than  in  isothermal  forging.  Clearly,  the  effect  of  tempera¬ 
ture,  possibly  combined  with  the  strain-rate  effect,  caused  the  metal  flow  to  differ 
under  the  two  forging  conditions.  Zienkiewicz  et.  al.  [33]  studied  and  performed 
a  coupled  thermal  analysis  for  steady-state  extrusion  operations.  Tang  et.  al. 

[34]  analyzed  the  shell  nosing  problem  using  this  approach.  Coupled  analysis 
using  the  updated  Lagrangian  approach  also  has  been  reported  m  earlier  work 

[35] . 

In  this  work,  the  finite  element  method  has  been  chosen  as  the  primary 
numerical  analysis  tool.  The  nonlinear  rigid  viscoplastic  finite  element  program 
ALPID  has  been  used  for  simulation  and  analysis  purposes  because  it  has  the 
following  capabilities: 

1.  Obtaining  detailed  solutions  of  mechanics  in  a  deforming  body  with  sufficient 
accuracy  for  practical  purposes.  The  solution  contains  velocities,  strains, 
stresses,  strain-rates,  temperatures,  contact  pressure  distributions,  die  load 
and  die  temperatures. 

2.  Handling  arbitrary  boundary  conditions. 
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3.  Including  the  friction  effect  at  the  die- workpiece  interface  in  both  deformation 
and  thermal  calculations. 

4.  Graphics  display  for  post-processing. 

5.  Analyzing  a  large  variety  of  problems  by  simply  changing  the  input  data. 

1.4  Design  of  optimal  process  parameters 

Compared  to  traditional  design  methods  in  the  metal  forming  field,  numer¬ 
ical  analysis  and  design  techniques  are  less  time  consuming  and  expensive.  In 
the  past,  considerable  research  has  been  conducted  in  designing  and  selecting 
optimum  process  parameters  using  numerical  approaches. 

Boer  et.  al.  [36]  developed  a  process  model  based  on  the  slab  method  to  calcu¬ 
late  stress  and  strain  distributions  in  the  deforming  workpiece.  Thermal  analysis 
was  performed  using  finite  element  methods  and  a  simplified  non-dimensional 
analysis.  Using  the  above  strategy,  and  minimizing  the  stress  ratio  parameter, 
optimal  ram  velocity  profiles  and  initial  die  temperature  were  obtained  for  both 
isothermal  and  hot  die  forgings  (using  NIM80A  material).  Lanka  and  Grand- 
hi  [37]  developed  the  Conformal  Mapping  Method  to  design  intermediate  die 
shapes  for  2-D  and  3-D  forgings.  This  is  a  geometry  mapping  technique  wherein 
the  staging  criteria  are  identified  based  on  the  stress  ratio  parameter.  Malas  [38] 
developed  an  approach  for  process  parameter  design  using  a  linear  relationship 
between  the  ram  velocity  and  strain-rate  as  an  approximation  to  maintain  the 


billet  variables  within  stable  processing  regions.  Hong  [39]  designed  a  control 
scheme  based  on  a  finite  element  analysis  model.  This  model  utilizes  the  local 
thickness  of  the  deforming  blank  as  a  measured  variable  to  track  a  desired  tra¬ 
jectory.  Kobayashi  et.  al.  [27],  and  Park  et.  al.  [40]  developed  the  “backward 
tracing”  technique  which  is  based  on  the  reversal  of  material  flow  during  simu¬ 
lation.  Using  this  method  preforms  were  designed  for  shell  nosing,  plane-strain 
rolling,  and  axisymmetric  forging  problems.  Han  et.  al.  [41]  combined  this  idea 
with  a  numerical  optimization  approach,  and  designed  optimal  intermediate  die 
shapes  for  isothermal  forging  processes  by  minimizing  the  strain-rate  variance  in 
the  deforming  workpiece.  This  technique  is  called  the  Backward  Deformation 
Optimization  Method”,  and  includes  sensitivity  analysis,  besides  introducing  a 
criterion  for  nodal  separation  from  the  surface  of  the  die  during  backward  de¬ 
formation  simulation.  Grandhi  et.  al.  [42]  then  introduced  an  optimal  control 
design  algorithm  into  the  process  parameter  design  procedure.  The  metal  form¬ 
ing  process  was  modeled  and  condensed  into  the  state  space  form,  and  a  suitable 
optimal  control  algorithm  was  used  in  designing  the  process  parameters.  Opti¬ 
mum  ram  velocity  schedules  for  maintaining  specified  strain-rates  in  the  billet 
were  generated  by  this  approach  for  an  isothermal  disk  forging. 

During  conventional  hot  die  forging,  there  is  a  complex  thermal  interplay 
between  the  workpiece,  die(s)  and  the  atmosphere.  There  is  heat  generation  due 
to  the  dissipation  of  deformation  energy,  and  friction.  At  the  same  time  there 
is  heat  loss  from  the  system  due  to  conduction  between  the  billet  and  die.  In 
addition,  there  is  also  heat  loss  due  to  convection  and  radiation  between  the  billet 
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and  the  environment.  As  a  result,  the  processed  material  may  experience  severe 
temperature  gradients  along  the  interface  of  the  die  and  billet,  which  could  lead 
to  the  formation  of  surface  defects.  Thermal  disparities  and  temperature  changes 
in  the  workpiece  may  also  induce  phase  transformations  and  changes  in  the  grain 
structure.  These  changes,  in  turn,  affect  the  flow  stress  and  metal  flow  as  well 
as  other  process  variables.  Furthermore,  severe  temperature  gradients  result  in 
large  thermal  stresses  leading  to  material  failure.  Quenching  is  an  example  where 
thermally  induced  stresses  can  cause  warping  and  cracking  of  the  finished  parts. 
In  view  of  these  effects,  this  research  extended  the  isothermal  study  in  reference 
42  to  handle  nonisothermal  situations  like  hot  die  forging. 

The  objective  of  this  work  is  to  build  a  numerical  model  representing  the 
metal  forming  system,  and  design  optimal  process  parameters  (ram  velocity  and 
initial  die  temperature)  that  satisfy  the  following  requirements: 

1.  Maintain  the  strain-rate  at  a  certain  value  at  a  given  location  (element)  in 
the  billet. 


2.  Maintain  the  temperature  at  the  critical  spot  (node)  above  some  value. 

3.  Reduce  the  temperature  range  in  the  billet. 

4.  Reduce  the  temperature  gradient  at  the  interface  of  die  and  billet. 

5.  Force  the  strain-rate  and  temperature  to  follow  desired  trajectories. 


* 


* 
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The  above  requirements  have  significant  practical  meaning,  as  explained  in 
the  later  sections  of  this  report.  If  these  requirements  can  be  satisfied,  maintain¬ 
ing  the  process  parameters  within  favorable  processing  regions  becomes  easy. 

The  coupled  deformation  and  thermal  analysis  code,  ALPID,  is  modified  to 
build  the  state  space  model  and  then  condensed  to  obtain  a  reduced  order  model. 
An  optimal  finite  time  controller  design  algorithm  has  been  integrated  with  the 
ALPID  code  to  calculate  the  required  (optimal)  process  parameters.  The  results 
show  that  the  above  methodology  is  quite  successful  in  achieving  the  required  ob¬ 
jectives.  But  developing  the  state  space  model  from  the  finite  element  equations 
and  solving  the  resulting  set  of  equations  is  a  non-trivial  task,  especially  while 
dealing  with  large  scale  systems.  While  simulating  realistic  manufacturing  pro¬ 
cesses,  it  is  necessary  to  use  large-scale  finite  element  models  with  a  large  number 
of  degrees  of  freedom.  In  such  instances,  the  corresponding  state  space  model 
also  has  a  large  number  of  states.  This  is  likely  to  cause  numerical  difficulties 
during  implementation,  besides  being  computationally  expensive  and  tedious.  In 
such  situations,  model  reduction  techniques  have  to  be  used  to  reduce  the  order 
of  the  system  before  designing  the  controllers.  The  following  section  gives  a  brief 
introduction  and  description  of  model  reduction  techniques. 


» 
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1.5  Reduced  order  models 


Balas  [43],  surveyed  different  model  reduction  methods  applied  to  the  con¬ 
trol  of  large  space  structures.  He  presented  a  mathematical  framework  to  survey 
the  general  topics  in  control  theory  (for  large  space  structures)  like  optimal  con¬ 
trol,  reduced-order  models,  spillover  compensation,  reduced  order  controllers, 
and  other  related  topics.  These  numerical  techniques  are  widely  used  in  aircraft 
and  space  structure  control  work  by  many  researchers.  Marshall  [44]  described 
in  detail  the  numerical  strategies  used  for  the  reduction  of  large  models  using  the 
modal  analysis  approach.  Akoi  [45,46]  applied  economic  aggregation  techniques 
to  solve  large  scale  control  problems.  Rao,  et.  al.  [47]  applied  the  balanced  trun¬ 
cation  and  modal  aggregation  methods  to  large  flexible  structures,  where  the  FE 
model  of  the  structure  in  the  state  space  form  is  used  to  design  robust  controller- 
s  for  structural  systems.  Safonov,  et.  al.  [48]  used  the  modal  truncation  and 
Hankel-norm  techniques  to  reduce  a  116-state  model  to  a  4-state  model,  while  de¬ 
signing  a  robust  multivariable  controller  for  suppressing  active  vibration  in  large 
space  structures.  Ben  Jaafar,  et.  al.  [49]  had  earlier  applied  model  reduction 
techniques  to  thermal  diffusion  problems.  In  that  work,  the  Eitelberg,  Marshal- 
1,  and  Aggregation  methods  were  applied  to  a  finite  element  model  describing 
heat  transmission  in  thermal  diffusion  problems.  Wanxie,  et.  al.  [50]  used  the 
relationship  between  the  generalized  variational  approach  and  LQ  control  theo¬ 
ry  to  successfully  reduce  the  size  of  an  eigenvalue  problem  by  half.  Chang  and 
Engblom  [51]  used  the  Rayleigh-Ritz  method  to  reduce  the  model  of  a  structure 
under  uniform  loading  for  non-linear  dynamic  response  predictions. 
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In  most  of  the  above  applications,  the  model  reduction  techniques  used  tend 
to  retain  the  dynamic  response  of  the  original  system.  For  the  current  analysis, 
the  objective  during  model  reduction  is  to  derive  a  reduced  order  system  that 
is  both  controllable  and  observable.  In  addition,  it  must  be  ensured  that  the 
reduced  model  retains  the  properties/characteristics  of  the  full  model.  Gregory 
[52]  constructed  a  reduced  model  for  lightly  damped  structures  using  balanced 
model  reduction  (BMR)  techniques.  He  reduced  a  full  state  model  having  114 
states  to  a  reduced  model  having  26  states,  and  concluded  that  the  BMR  tech¬ 
nique  (effectively)  retains  the  stability,  controllability  and  observability  of  the 
full  model.  Yae  and  Inman  [53],  also  concluded  that  the  BMR  technique  retains 
the  properties  of  the  full  state  model  for  multi-body  systems  and  can  be  used  for 
control  applications.  This  was  also  confirmed  by  Adams  et.  al.  [54],  who  used 
the  BMR  method  in  flight  control  applications. 

In  this  work,  the  following  three  broad  methods  of  model  reduction  were 
studied: 

a.  Aggregation  Method  [55] 

b.  Modal  analysis  methods  [56] 

c.  Balanced  model  reduction  [57] 

After  the  reduced  order  state  space  models  are  developed,  a  control  law  based 

on  the  differential  Riccatti  equation  is  adopted  to  design  optimal  die  velocities  for 
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the  metal  forming  process,  like  in  the  case  of  the  full  state  model.  The  detailed 
formulation  involved  in  deriving  reduced  order  models  is  explained  at  a  later 

stage  in  this  report. 

The  following  chapters  describe  in  detail  the  process  control  strategy  devel¬ 
oped  and  the  mathematical  formulations  involved  in  setting  up  and  solving  the 
optimal  control  problem.  Appendix  B  and  Appendix  C  provide  information  on 
the  program  developed  (COPP  -  Control  of  Optimal  Processing  Parameters)  for 
implementing  this  methodology. 


CHAPTER  2 


State  Space  Model  for  Deformation  Analysis 


2.1  Basis  for  the  finite  element  formulation 

The  finite  element  method  (FEM)  used  in  this  work  is  based  on  the  flow 
formulation.  This  approach  assumes  that  during  metal  forming  the  plastic  s- 
trains  usually  outweigh  elastic  strains  and  the  idealization  of  rigid-plastic  or 
rigid-viscoplastic  material  behavior  is  acceptable.  In  other  words,  phenomena 
associated  with  elasticity  are  completely  neglected  in  this  method  of  analysis. 

The  original  problem  associated  with  the  deformation  process  of  rigid  vis¬ 
coplastic  materials  is  a  boundary-value  problem  and  the  formulation  is  briefly 
described  below  [22,23,27]: 

Consider  a  body  (workpiece)  having  volume  V  and  a  boundary  surface  S. 
The  boundary  surface  S  may  be  divided  into  three  distinct  parts  given  by: 

S  =  SU  +  SF  +  SC  (2.1) 

where  Su  is  the  surface  with  prescribed  velocity  it,  Sp  is  the  portion  of  the  work- 
piece  with  traction  F  prescribed,  and  Sc  is  the  tool-workpiece  interface  surface 
where  the  frictional  stress  /  acts.  The  body  is  composed  of  a  rigid  plastic  mate¬ 
rial  which  obeys  the  Von  Mises  yield  criterion  and  its  associated  flow  rule.  The 
deformation  of  the  body  V  is  now  characterized  by  the  following  field  equations. 
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Equilibrium  conditions  for  deformation  analysis  (neglecting  body  forces)  are 
given  by: 

=  o  (2.2) 

where  <x;;-  is  the  ijth  component  of  the  Cauchy  stress  tensor.  In  the  above  equa¬ 
tion,  a  recurring  letter  suffix  indicates  the  sum,  and  the  comma  denotes  partial 
differentiation. 


Compatibilty  conditions  are  given  by  the  following  strain  rate- velocity  rela¬ 


tion: 


iij  =  2  (“W  +  ui. •) 


(2.3) 


where,  tij  is  the  ij(*fe-component  of  the  strain-rate  tensor,  u;  is  a  velocity  compo¬ 
nent,  and  the  comma  represents  differentiation. 

The  constitutive  equation  giving  the  stress-strain  rate  relationship  is  [27]: 

.  2  cr . 


aii  “3  iZii 


(2.4) 


where  cr'ij  is  a  component  of  the  deviatoric  stress,  a  and  e  are  the  effective 
stress  and  effective  strain  rates,  respectively,  defined  by  a  =  o’ij'  and 


t  =  \J\kijkij.  The  flow  stress,  in  general,  is  a  function  of  total  strain,  strain-rate, 
temperature,  and  microstructure  S,  and  may  be  represented  as: 


v  =  f(T,e,k,S) 


(2.5) 


The  boundary  condition  for  this  system  representing  equilibrium  of  stresses  is 
given  by: 

crijTii  —  Fj,  on  Sf  (2-6) 
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where  n{  is  the  ith-component  of  the  unit  normal  vector  to  the  body  surface 
and  Fj  is  the  jth-component  of  the  prescribed  surface  traction.  The  velocity 
boundary  condition  is  given  by: 

u{  =  i 2{,  on  Su  (2-7) 

where  Ui  is  the  ith-component  of  prescribed  velocity.  Solutions  to  this  problem 
are  the  stress  and  velocity  distributions  that  satisfy  the  governing  equations 
and  the  boundary  conditions.  The  governing  equations  include  the  equilibrium 
equations,  the  yield  criterion,  and  the  compatibility  conditions  derived  from  the 
flow  rule.  Since  it  is  difficult  to  obtain  a  complete  solution  that  satisfies  all 
the  field  equations,  various  approximate  methods  are  used  to  solve  the  above 
problem,  one  of  them  being  the  finite  element  method.  The  basic  principles  and 
concepts  involved  in  the  finite  element  method  are  the  variational  principle  and 
discretization. 

The  finite  element  governing  equation  for  metal  forming  analysis  may  be 
derived  from  the  potential  energy  functional  [27], 

7t  =  f  aidV  -  [  FmdS  (2.8) 

Jv  JSf 

where  the  first  term  represents  the  rate  of  total  plastic  work  done  and  the  sec¬ 
ond  term  represents  the  rate  of  work  expended  due  to  surface  traction.  The 
variational  principle  requires  that  among  admissible  velocities  that  satisfy  the 
conditions  of  compatibility  and  incompressibility  (as  well  as  the  velocity  bound¬ 
ary  conditions)  the  actual  solution  gives  the  above  functional  a  stationary  value. 
In  other  words,  the  solution  of  the  original  boundary  value  problem  is  obtained 
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from  the  solution  of  the  dual  variational  problem,  where  the  first  order  variation 
of  the  above  functional  vanishes,  i.e., 

—  [  aStdV  -  /  FiSmdS  =  0  (2.9) 

JV  JSf 

* 

Metal  deformation  occurs  at  a  constant  volume.  Therefore,  the  admissible 
velocity  U{  needs  to  satisfy  the  incompressibility  constraint.  This  is  embedded 
into  the  variational  principle  functional  (Eq.  (2.9))  by  introducing  the  penalty 
constant  K  (which  is  a  very  large  positive  quantity)  into  the  equation  as, 

8tt  =  [  <j8edV  -f  K  [  e'vSe'vdV  —  [  Fi8u{dS  =  0  (2.10) 

JV  Jv  J  Sp 

where  by  =  e«  is  the  volumetric  strain  rate.  The  development  of  the  finite 
element  equations  from  the  above  equation  is  described  in  the  next  section. 

2.2  Finite  element  equations 

The  discretization  of  the  variational  principle  functional  decribed  above  is 
done  using  the  standard  procedure  of  the  finite  element  method.  The  primary 
unknown  for  the  solution  of  Eq.  (2.10)  (which  represents  a  quasi-static  plastic 
deformation  process)  is  the  velocity  field  associated  with  it.  This  velocity  field, 
u ,  is  approximated  by  shape  functions  in  terms  of  nodal  point  velocity  values  as, 

u  =  Npv  (2-11)  P 

where  Nd  is  the  shape  function  matrix  for  deformation  analysis,  and  v  is  the 

m 

nodal  velocity  vector  for  the  element  under  consideration.  For  two-dimensional, 
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four-noded,  quadrilateral  elements  having  two  degrees  of  freedom  per  node,  v  is 
a  8  x  1  vector  given  by, 

VT  =  [vi  V2  V3  .  ..Us] 

Eq.  (2.10)  may  now  be  expressed  in  terms  of  the  nodal  point  velocities  v  and 
their  variations  6v.  Using  the  variational  principle,  from  the  arbitrariness  of  Svj, 
a  set  of  nonlinear  algebraic  equations  (stiffness  equations)  are  obtained  as  given 
below  [27]: 


where  ( j )  represents  the  quantity  at  the  jth  element.  The  capital-letter  suffix, 
/,  refers  to  the  nodal  point  number.  The  above  equation  is  actually  determined 
by  first  obtaining  the  elemental  equations  and  then  assembling  them  under  ap¬ 
propriate  constraints  to  obtain  a  set  of  globed  stiffness  equations  in  the  standard 
finite  element  form, 

nK"V  =  (2.13) 

where  K  is  the  global  stiffness  matrix,  F  is  the  load  vector  (described  in  detail  in 
the  next  section),  V  is  the  global  velocity  vector,  and  n  is  the  counter  representing 
the  iteration  number  at  every  time  step.  This  problem  is  similar  to  a  standard 
finite  element  analysis  problem,  but  its  special  feature  is  that  the  geometry  of  the 
boundary,  and  hence  the  boundary  condition,  keeps  changing  as  the  die  stroke 
increases.  K  and  F  are  both  functions  of  V,  resulting  in  a  highly  nonlinear 
system  of  equations.  Therefore,  the  analysis  path  of  the  forging  process  is  traced 
in  an  incremental  manner  by  considering  a  series  of  discrete  equilibrium  states, 
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each  corresponding  to  a  specific  value  of  time.  The  discrete  equilibrium  state 
is  assumed  to  correspond  to  the  time  value  t,  and  the  subsequent  or  (n  +  1) 
equilibrium  state,  is  assumed  to  correspond  to  time  ( t  +  At).  The  geometry 
configuration  at  time  ( t  +  At)  is  updated  from  time  (t).  This  is  known  a  priori 
from  the  solution  of  the  previous  time  step  (t).  The  velocity  distribution  for 
configuration  (l t  +  At)  is  unknown  and  calculated  in  an  implicit  iterative  manner 
until  a  converged  solution  is  obtained. 


The  solution  to  Eq.  (2.12)  or  Eq.  (2.13)  is  generally  obtained  iteratively 
using  the  Newton-Raphson  method.  The  method  consists  of  linearization  and 
application  of  convergence  criteria  to  obtain  the  final  solution.  If  the  converged 
point  is  v  =  vq,  a  Taylor  series  expansion  can  be  made  about  Vq  as  follows: 


dir 


dvj 


+ 


v=v0 


d2 


7T 


Avj  =  0 


(2.14) 


V=v0 


[dvjdvj 

where  Avj  is  the  first  order  correction  of  the  velocity  v.  The  linearized  system 
of  equations  may  be  written  as: 


KAv  =  f 


(2.15) 


where  K  is  again  the  material  and  process  dependent  stiffness  matrix,  and  f 
is  the  residual  of  the  nodal  point  force  vector.  This  equation  may  further  be 
rearranged  and  expanded  as, 

Ksv  +  KtAv  =  F  (2.16) 

where  K s  and  K t  are  called  the  secant  stiffness  and  tangent  (or  gradient)  stiffness 
matrices,  respectively  (Fig.  1).  K5  and  K*  are  calculated  based  on  the  current 


» 


* 


22 


solution  of  Eq.  (2.13)  and  treated  as  constant  matrices  for  time  t  to  (t  +  At).  Here 
again,  F  is  the  finite  element  load  vector.  Once  the  solution  of  Eq.  (2.14)  for  the 
velocity  correction  term  (Av)  is  obtained,  the  current  velocity  Vo  is  updated  as 
(Vo+oAv),  where  a  is  a  constant  between  0  and  1,  and  is  called  the  deceleration 
coefficient.  Iteration  is  continued  until  the  velocity  correction  terms  become 
negligible  and  the  convergence  criteria  are  met. 

To  calculate  the  first  and  second  derivatives  in  Eq.  (2.14),  the  effective  strain 
rate  (I)  and  volumetric  strain  rate  (e^)  must  be  expressed  in  terms  of  the  nodal 
velocities.  This  is  done  by  means  of  the  strain-rate  matrix  Bs  defined  by  [27], 

e  =  Bsv  (2-17) 

where  e  is  the  strain-rate  vector  which  contains  the  normal  strain-rate  and  engi¬ 
neering  shear  strain-rate  components.  The  effective  strain-rate,  in  turn,  can  be 
expressed  in  terms  of  the  strain-rate  vector  as  [27]: 

(I)2  =  eTD£ 

=  vTBjDBsv  (2-18) 

=  vTPv 

where  P  =  BjDBs.  D  is  a  diagonal  matrix  that  has  |  and  |  as  its  compo¬ 
nents,  corresponding  to  the  normal  and  shear  strain-rates,  respectively.  Also, 
the  volumetric  strain-rate  is  defined  as: 

tv  =  Cm  (2.19) 
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where  Ci  =  Bu  +  B^i  +  B 3;,  and  Bij  is  an  element  of  the  strain-rate  matrix  Bs. 
Using  the  above  formulae,  the  first  and  second  order  derivatives  in  Eq.  (2.14) 
may  be  expressed  as  [27], 


dir 

dv 


-  =  [  lpijVjdV  +  [  KCjVjCidV  -  [  FjNjidS  (2.20) 

;  Jv  e  JV  JSf 


and, 

927T 

dv{Vj 


=  jv  | PijiV  +  f( -  p)  IPwmP^dV  +  jv  KCjCiiV  (2.21) 


The  detailed  derivation  of  the  above  equations  and  explanation  of  the  terms 
involved  may  be  obtained  from  reference  27. 


2.3  Frictional  force  as  the  finite  element  load  vector 


The  frictional  condition  at  the  die-material  interface  greatly  influences  metal 
flow  during  most  forming  processes  and  merits  proper  analysis  and  treatment. 
In  the  current  procedure,  friction  stress  //,  which  appears  in  the  load  vector  of 
the  finite  element  governing  equation,  is  given  by, 


(2.22) 


where  m  is  the  friction  factor,  k  is  the  shear  strength  of  the  material,  and  uT  is 
the  sliding  velocity  at  the  die  workpiece  interface  (Sc)- 

During  finite  element  analysis,  a  coordinate  transformation  is  made  to  trans¬ 
fer  the  die-workpiece  contact  nodal  velocities  onto  a  local  coordinate  system 
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which  has  one  axis  normal  to  the  die  surface  and  the  other  axis  along  the  die 
surface  (Fig.  2).  A  shape  function  is  used  to  interpolate  ur  [58]  as, 


Ur  =  ^  (1  -  r)  vr(j )  +  ^  (1  +  r)  ur(/)  (2.23) 

where  vT^  and  vr ^  are  the  relative  sliding  velocities  at  nodes  j  and  1,  respectively. 
At  node  j,  the  sliding  velocity  can  be  expressed  in  the  local  coordinate  system 
as, 

vT(j)  =  vT(j)  ~  vdT(j)  (2-24) 


where  is  the  tangential  nodal  velocity  along  the  die  surface,  and  ^dT(j)  is  the 
tangential  component  of  the  die  velocity  at  boundary  node  j.  Then,  Eq.  (2.22) 
becomes: 


*  =  -  bmkNTf^\iS 


(2.25) 


~L 


Ur 


L 


mkrf^dS 


uT 


where  Ny  is  the  shape  function,  vy  and  v^y  are  given  below: 


y 

vT 


[VT(;)  vT(o]  VJj  -  [VdT(j)  vdT(o] 


Eq.  (2.25)  gives  the  friction  force  along  the  die  surface  for  the  segment  between 
node  j  and  node  1.  The  first  term  of  Eq.  (2.25)  is  added  to  the  stiffness  matrix 
K3  because  it  contains  the  nodal  velocity.  The  second  term  stays  on  the  right 
side  of  the  finite  element  governing  equation,  and  appears  as  the  load  vector  in 
Eq.  (2.16).  During  the  design  phase,  this  term  is  transformed  into  the  control 
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Workpiece  Surface 


n  s  =  Inward  unit  normal  to  the  die  continuum  at  the  contact 
node  i 

n  =  Unit  vector  tangent  to  the  die  continuum  at  the  contact 
node  i 

(  n  s  5  n  r  )  =  Local  coordinates 
(  X  ,  y  )  =  Global  coordinates 


Figure  2:  Contact  boundary  conditions 
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input  as  it  contains  die  velocity,  which  is  one  of  the  major  design  parameters  in 
this  work. 

2.4  Condensation  of  the  system  during  deformation  analysis 

In  general,  while  simulating  metal  forming  processes,  large-scale  finite  ele¬ 
ment  models  are  required  to  capture  in  detail  the  thermo-mechanical  behavior 
of  the  deforming  material.  In  this  work  four-noded  isoparametric  quadrilateral 
elements  have  been  used  for  discretization  of  the  workpiece  and  the  die.  Each 
of  these  nodes  has  two  degrees  of  freedom.  Because  each  nodal  degree  of  free¬ 
dom  corresponds  to  one  algebraic  equation,  the  number  of  equations  to  be  solved 
during  analysis  would  be  very  large  if  many  elements  are  used  to  describe  the 
problem.  This  would  result  in  a  large  size  state  space  model,  adding  to  the  dif¬ 
ficulty  and  computational  expense  in  integrating  FEM  with  an  optimal  control 
algorithm.  In  addition,  most  control  design  algorithms  are  efficient  and  effec¬ 
tive  only  while  handling  50  or  less  state  variables.  Beyond  this  limit,  numerical 
difficulty  may  be  encountered  during  the  computation  and  implementation  pro¬ 
cess.  There  is  thus  a  need  for  reducing  the  nodal  degrees  of  freedom  of  the  metal 
forming  system  represented  by  Eq.  (2.16). 

A  systematic  condensation  procedure  has  been  developed  wherein  some  of 
the  important  degrees  of  freedom  of  the  system,  including  the  nodal  degrees  of 
freedom  for  the  ‘element  of  interest’,  are  retained,  and  those  for  the  rest  of  the 
elements  are  condensed  out  of  the  system.  ‘Element  of  interest  actually  refers 
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to  the  element (s)  constituting  critical  regions  in  the  billet  where  it  is  desired  to 
achieve  and  maintain  specific  processing  conditions. 

To  condense  Eq.  (2.16),  the  velocity  field  V  has  been  divided  into  the  fol¬ 
lowing  five  zones  (Fig.  3): 

VT:  Tangential  nodal  velocities  along  the  die  surface. 

Vjv’:  Nodal  velocities  with  direction  normal  to  the  die  surface. 

Vc:  Fixed  nodal  degrees  of  freedom  due  to  symmetric  boundary  condition. 

Vj:  Nodal  velocities  of  the  element(s)  of  interest. 

V#:  Nodal  velocities  of  the  other  elements. 


Eq.  (2.16)  can  now  be  rewritten  based  on  the  above  definitions  as, 


KsU  Ksi2  Ka\z  Ks  14  Ks  15 

-VT- 

K32\  Ks22  K3 23  K3  24  Ks  25 

VN 

K3  31  K3 32  Ks33  Ks34  -^>35 

Vi 

Kail  K3  42  Ks  43  K3  44  K3a5 

Vb 

_Ks  51  KsS2  Ks53  Ks54  Ks  55. 

LVcJ 

Km  Ktl2  KtlZ  Kti4  Ktl5 

'  AVx 

Ft 

Kt2l  Kt22  Kt23  Kt2i  Kt25 

AVn 

Fn 

Kfi\  Ki32  Ki33  Kt34  Kt35 

AVX 

= 

0 

Kt4i  Kt42  K%43  Km  Kt45 

AVb 

0 

.  KtSl  Kf52  Kt53  KtS4  KtSS . 

AVC. 

.  0  . 

where  K„;  and  Ktij  are  submatrices  of  K,  and  Kt,  and  are  determined  based 
on  their  location  in  the  billet.  FT  is  the  load  vector  associated  with  the  friction 
force,  and  FN  is  the  load  vector  related  to  the  normal  force  at  the  die/billet 

interface. 


* 
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Bottom  Die 


•  =  Nodes  at  the  die  contact  boundary  (  V  t  and  Vn  ) 

■  =  Nodes  at  the  symmetric  boundary  (Vc) 

A  =  Nodes  of  the  element  of  interest  (  Vi  ) 

# 

o  =  Nodes  of  the  rest  elements  (  V  B  ) 

Figure  3:  System  condensation 


30 


Assuming  die  velocity  is  represented  by  the  tangential  and  normal  die 
velocities  at  the  contact  boundary  are: 


1 

Normal 

Vdjy  =  T  yVd 

(2.27) 

Tangential 

V  dT  =  TzVd 

(2.28) 

i 


where  Tz  and  Ts  are  the  transformation  vectors  which  contain  the  direction 
cosines  of  the  local  coordinates  at  every  die  contacting  node.  Referring  to  Fig. 
2,  at  the  contact  node  i,  the  transformation  vectors  are  given  by, 


Txi  —  SZTl  OC{ 

(2.29) 

Tyi  =  cos  ai 

(2.30) 

In  the  nodal  velocity  discretization  scheme  described  above,  VN  is  always  equal 
t°  VdN,  the  normal  component  of  the  die  velocity.  Therefore,  it  is  reasonable 
to  set  AVjv  to  zero.  At  surfaces  having  symmetric  boundary  conditions,  the 
velocity  components  are  zero  in  the  direction  of  symmetry.  So,  we  also  have  Vc 
and  AVc  equal  to  0.  These  terms  can  therefore  be  dropped  from  the  system 
equations.  We  are  then  left  with  the  following  boundary  conditions: 

Vc  =  0  AVC  =  0 
VN  =  VdN  =  T  yVd  AVy  =  0 

By  applying  these  boundary  conditions  in  Eq.  (2.26)  and  deleting  the  corre¬ 
sponding  rows  and  columns,  we  have, 


Ksu  -Ksi3  KsU 
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Ktu  Km  Ktu 

AVt 

Ks31  Ks33  KsZA 

Vi 

+ 

Kfi\  Kt33  Kt34 

AVi 

Ks41  Ks43  KsU 

[VB. 

_  Km  Km  Kt44 . 

AVb. 
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(2.31) 


Ft  —  Ksi2Vn 
=  -K,  32VN 

-K,  42VN 

On  the  right  side  of  Eq.  (2.31),  Ft  is  obtained  from  the  friction  force  calculation 
of  Eq.  (2.25)  as, 


=  [  mkNTf^dS  Vd  (2’32) 

iSc  1  |uT| 

=  h3Vd 

where  V«,t  =  TxVd ,  and  ha  =  /Sc  mfcNj^dS.  Evaluating  the  above  friction 
term  along  the  whole  interface,  we  have, 

Ft  =  K3Vd  (2.33) 


The  load  vector  may  now  be  written  as, 

Ft  -  K8i2Vn  H3  —  Ks\2Ty 
—KsmVii  —  —KsZiTy  Vd  (2.34) 

Ks 42Vn  L  k«*Tv  - 

To  solve  and  condense  Vjg  from  Eq.  (2.31),  a  forward  difference  approach  is  used 
to  calculate  A Vb,  where, 


AVb  =  Vb  —  ~V  B  (2.35) 


Here,  \B  is  the  nodal  velocity  at  time  ( t  +  A t)  and  Vb  is  the  nodal  velocity 
at  time  t.  Vg  is  a  known  vector  and  treated  as  a  constant  vector  in  the  time 
interval  (t,t  +  At).  By  substituting  Eq.  (2.35)  into  the  3rd  row  of  Eq.  (2.31),  we 
have, 

Vb  =  [KS44  -f  Kt44] 
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[— Ks4iVt  —  Ks43 Vj  -  Kt4lAVT  —  Kt43AVi  —  Ks42 V n  +  Kt44Vfi] 

Substituting  this  solution  into  the  first  and  second  rows  of  Eq.  (2.31)  gives  the 
following  reduced  order  equation: 

' K3 ii  K3i2]\Vt]  \Ktu  Ktu]  [AVt]  _  [Gil  v  [Wil  (t) 

Ka 21  Ks22\  [vI\  +  [Km  Kt22\  [AVj\  -  [g2\  Vd^  [W2\  [Z-60) 

where  the  submatrices  and  vectors  are: 

K3n  =  Ksll  —  (Ksi4  +  KtX4)  (KS44  +  Kt44)  KS41 

Ks12  =  KS13  —  (KS14  +  Ktl4)  (K844  +  Kt44)  1KS43 
K321  =  Ks3i  —  (KS34  +  Kt34)  (K844  +  Kt44)  1KS4i 
K322  =  KS33  —  (KS34  +  Kt34)  (KS44  +  Kt44)  1KS43 
Kill  =  Ktll  —  (Ksi4  +  Kt14)  (KS44  +  Kt44)  XKt4i 
K«2  =  KtlS  —  (KS14  -f  Km)  (KS44  +  Kt44)  XKt43 

Km  =  Ktsi  —  (Ksi4  +  Kti4)  (KS44  +  Kt44)  XKt4i 

Kf22  =  Kt33  —  (KS14  +  Ktl4)  (KS44  +  K444)  XKt43 

Gi  =  Hs  —  Ksi2Ty  +  (K8i4  +  Km)  (KS44  +  Kt44)  1Ks42Ty 

G2  =  — K832Ty  +  (K834  +  Kt34)  (K844  +  Kt44)  KS42Ty 

Wi  =  -  (KS14  +  Ktu)  (K844  +  Kt44)  XKt44Vs  +  Ktl^Vs 
W2  =  -  (Ks34  +  Kt34)(K844  +  Kt44)  XKt44VB  +  Kt34Vfl 

In  Eq.  (2.36),  only  Vy  and  Vj  are  retained.  These  are  the  nodal  degrees  of 
freedom  needed  in  the  process  parameter  design.  The  rest  of  the  nodal  degrees 
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of  freedom  are  either  eliminated  or  condensed.  In  this  way,  the  order  of  the  state 
space  system  is  greatly  reduced. 


jt 


2.5  State  space  model  for  deformation  analysis 


In  Eq.  (2.36)  if  Ks,  K*,  G  and  W  are  defined  as  follows, 

T>  _  fKsii  KS12\ 
s  "  VKszi  ks22; 

it  —  ?tl2^ 

'  “  VKtai  KW 


and  if, 

VST  =  [vtt  V/]  A VST  =  [AVTT  AV7t] 

then  equation  (2.36)  may  be  written  as, 


KaVs  +  Kt  AVS  =  GVd  +  W 


or, 


Kt  AV5  =  -KaVs  +  GVd  +  W 


By  multiplying  and  dividing  the  left  side  of  the  above  equation  by  At,  we  have, 


[AiK,]  =  -K ,VS  +  GVd  +  W 

At  0 

Now,  by  applying  the  limit,  A t  — >  0,  we  have, 

[AtKt]  ^  =  -KsV5  +  GVd  +  W  (2.37) 
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where  A t  is  the  time  increment  from  one  simulation  step  to  the  next.  The  above 
equation  is  in  the  standard  state  space  format.  By  defining  the  plant  matrix  as, 

AD  =  -[AtKf]_1  [Ka] 

the  input  matrix  as, 

B  d  —  [AtKtj  [G] 

and  the  constant  perturbation  vector  (due  to  condensation)  as, 

Wd  =  [AticJ-1  [w] 

the  condensed  state  space  model  for  deformation  analysis  (Eq.  (2.37))  may  be 
rewritten  as: 

^  =  Ad  Vs  +  BDVd  +  Wd  (2-38) 

dt 

where  Vs  is  the  state  vector  containing  the  nodal  velocities  associated  with  the 
element  of  interest  and  the  nodal  velocities  of  the  die-contacting  boundary  nodes. 
The  plant  matrix,  input  matrix,  and  perturbation  vector  have  been  described 
above,  and  Vd  (die  velocity)  is  the  input  (scalar)  to  the  system  represented  in  the 
above  equation. 
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CHAPTER  3 


State  Space  Model  for  Thermal  Analysis 


3.1  Finite  element  equations  for  heat  transfer 


* 


The  consideration  of  temperature  effects  in  metal  forming  is  very  important 
because  thermal  effects  accompanying  the  deformation  process  strongly  influence 
the  mechanical  and  material  properties  of  the  final  product.  In  addition  to  defor¬ 
mation  analysis,  nonisothermal  metal  forming  simulation  needs  a  comprehensive 
thermal  analysis  procedure  to  determine  the  temperature  distribution  in  the  bil¬ 
let  and  die  domains.  The  thermal  analysis  involves  several  separate  bodies  like 
the  die,  workpiece,  and  lubricant,  and  takes  into  account  the  thermal  interplay 
between  each  of  them.  The  die  and  billet  are  discretized  (Fig.  4),  and  finite 
element  analysis  for  each  of  these  bodies  is  conducted  separately.  Heat  transfer 
between  the  distinct  bodies  through  the  region  of  contact  is  modeled  by  enforc¬ 
ing  consistent  heat  transfer  boundary  conditions.  Generally,  for  nonisothermal 
forming  processes,  a  coupled  thermo- viscoplastic  analysis  is  carried  out  wherein 
it  is  necessary  to  simultaneously  solve  the  material  flow  problem  (for  a  given 
temperature  distribution)  and  the  heat  transfer  equations.  This  section  gives  an 
overview  of  the  thermal  analysis  procedure  for  metal  forming  simulation  using 
the  finite  element  method  [27,34]. 


-4 
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Consider  a  body  in  thermal  equilibrium  under  a  specific  set  of  prescribed 
thermal  boundary  conditions  (Fig.  5).  The  energy  balance  equation  for  heat 
transfer  analysis  is  given  by, 

V  •  q  t  pc~~  —  0  (3.1) 

at 

where  p  and  c  are  the  mass  density  and  specific  heat,  respectively;  f  is  the  heat 
generation  rate,  T  is  the  current  temperature,  and  t  refers  to  time.  Also,  q  is 
the  heat  flux,  and  can  be  written  using  Fourier’s  law  of  heat  transfer  as, 

q  =  — kc  •  VT  (3.2) 

where  kc  denotes  thermal  conductivity,  and  VT  is  the  spatial  gradient  of  tem¬ 
perature.  In  Eq.  (3.1),  V  •  q  refers  to  the  divergence  of  heat  flux,  which  is  given 
by,  V  •  q  =  where  Xm  is  the  m-coordinate  value  for  a  generic  point  P  in 

the  body  domain.  If  the  heat  generation  in  the  deforming  body  is  assumed  to  be 
due  to  plastic  deformation  only,  then, 

T  =  K<TijSij  (3-3) 

where  the  heat  generation  efficiency  k  represents  the  fraction  of  mechanical  en¬ 
ergy  transformed  into  heat  (which  is  generally  assumed  to  be  0.9).  In  this  work, 
two  types  of  boundary  conditions  are  considered  over  the  body  surface  (Fig.  5). 
On  surface  St,  the  temperature  is  prescribed  as: 


Interface  Hat  Transfer  Between 
the  Specirraen  and  Die 


i 


Convection  and 
Radiation  Boundary 
Conditions 


Prescribed  Temperature 


Define: 

Sc  =  Portion  of  Sq  where  convection  and  radiation  boundary 
conditions  are  specified. 

S  t  =  Portion  of  Sq  where  interface  heat  transfer  conditions  are  applied. 
S  q  =  S  c  +  S  i 


1 


Figure  5:  Thermal  boundary  conditions 
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On  surface  Sq,  the  heat  flux  is  specified  as: 

n  ■  q  +  q  =  0  on  Sq  —  S]  -f  Se  (3*5) 

where  q  is  the  heat  flux  due  to  conduction,  and  n  is  a  unit  vector  acting  normal 
to  the  body  surface  in  an  outwards  direction,  q  represents  heat  flux  due  to 
convection  (and/or  radiation)  boundary  conditions  and  heat  generated  due  to 
friction  in  the  die-workpiece  contact  area.  Sj  is  the  die-workpiece  interface,  and 
se  represents  the  surface  exposed  to  the  environment. 

Using  the  Galerkin  weighted  residual  approach,  with  ST  (i.  e.,  virtual  temper¬ 
ature  increment)  as  the  weighting  function,  a  weak  form  of  Eq.  (3.1)  is  obtained 
as  [27], 

Jv  |v  ■  (kc  •  VT)  +  r  —  Pc~^  6TdV  +  Js(q  +  n-q)6TdS  =  0  (3.6) 

where  S  is  the  total  surface  area  of  the  body  ( S  =  St  -f  Sq),  and  V  is  the 
current  volume  of  the  body.  Solutions  to  problems  of  this  nature  require  that 
the  temperature  field  satisfies  the  prescribed  boundary  conditions  and  Eq.  (3.6) 
for  an  arbitrary  perturbation  St. 

To  implement  the  finite  element  procedure,  the  temperature  field  in  Eq.  (3.6) 
is  discretized  and  approximated  as, 

T  =  NjT  (3.7) 

For  a  quadrilateral  element,  the  shape  function  vector,  N,  and  the  vector  of  nodal 
temperatures,  T,  are  given  by, 

N T  =  [qi  72  73  74]  TT  =  [Ti  r2  T3  T4] 
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where  $  is  the  shape  function,  and  T;  is  the  nodal  temperature.  By  substituting 
Eq.  (3-7)  into  Eq.  (3.6)  and  expanding,  the  finite  element  governing  equation 

for  heat  transfer  analysis  is  obtained,  as  shown  below: 

CtT  +  KCT  =  Qt  (3-8) 

where  CT  is  the  heat  capacity  matrix,  Kc  is  the  heat  conduction  matrix,  Qt  is 
the  heat  flux  vector,  and  T  is  the  vector  of  nodal  point  temperatures. 

The  transient  thermal  behavior  of  the  body  is  captured  in  an  incremental 
manner  by  considering  a  series  of  discrete  thermal  equilibrium  states,  each  corre¬ 
sponding  to  a  specific  value  of  time.  The  nth  discrete  equilibrium  state  is  assumed 
to  correspond  to  the  time  value  (t)  and  the  subsequent,  or  (n  +  1)*A  equilibrium 
state  is  assumed  to  correspond  to  time  (t  +  A t)  where  A t  is  the  (current)  time 
increment  value.  Therefore,  for  a  particular  CT  and  Kc,  Eq.  (3.8)  is  used  to 
obtain  the  temperature  solution  in  the  interval  (t,t  +  At).  Then,  CT  and  Kc  are 
updated,  and  the  procedure  is  repeated  for  the  next  time  step. 

In  Eq.  (3.8),  the  heat  flux  vector  Qt  has  several  components  and  is  expressed 
as: 

Qt  =  f  K,(ae )  N t^V  +  f  ere  (T f  —  Ts )  N jdS 
JV  JSe 

(3.9) 

+  [  h(Te-T,)NTiS+  [  h,ui(Td-T„)NTiS+  I '  qfNTdS 

JSe  JSI  1 

where  Se  refers  to  the  surface  where  convection  and  radiation  boundary  con¬ 
ditions  are  specified.  5/  is  the  surface  where  interface  heat  transfer  conditions 

apply. 
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The  first  term  on  the  right  side  of  Eq.  (3.9)  is  the  heat  generated  by  plastic 
deformation  within  the  deforming  body.  The  second  term  defines  the  heat  ra¬ 
diated  between  the  workpiece  and  environment,  where  cr  is  the  Stefan-Boltzman 
constant,  e  is  the  emissivity,  and  Te  and  Ts  are  environment  and  surface  tem¬ 
peratures,  respectively.  The  third  term  describes  the  heat  convected  from  the 
body  to  the  environment  with  heat  convection  coefficient  h.  The  fourth  term 
represents  the  heat  conducted  between  the  die/workpiece  through  their  inter¬ 
face.  Td  and  Tw  are  die  and  workpiece  temperatures,  respectively,  and  h[u j,  is  the 
heat  transfer  coefficient  for  the  lubricant.  The  last  term  is  the  contribution  of 
the  heat  generated  due  to  friction  along  the  die- workpiece  interface,  qf  being  the 
surface  heat  generation  rate  due  to  friction. 

3.2  State  space  model  for  thermal  analysis 

From  the  finite  element  governing  equations  for  thermal  analysis,  a  state 
space  model  representing  the  thermal  aspects  of  the  metal  forming  process  is 
constructed.  To  build  the  state  space  model  for  thermal  analysis,  Eq.  (3.8)  first 
must  be  expressed  in  terms  of  the  state  variables,  namely  the  nodal  velocities  and 
nodal  temperatures.  The  most  difficult  task  in  achieving  this  is  the  representation 
of  the  thermal  load  vector  of  Eq.  (3.8)  in  terms  of  the  state  variables.  This  is  done 
by  considering  one  term  at  a  time  from  Eq.  (3.9).  The  procedure  is  described 
below: 
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(1)  fvn(ai)NTdV 

Using  Eq.  (2.18)  to  describe  £,  we  have, 


J  k(<x£)N 'idV  =  N EdV 


L 


^NTvTP  dVv 

V  £ 


(3.10) 


=  [  “Nrvjp<n'v 

JV  £ 


=  KEV 

where  v  is  the  nodal  velocity  vector  for  the  element  under  consideration  (un¬ 
known),  K E  =  fv  and  v*  is  the  vector  of  nodal  velocities  (known 

quantity)  at  the  beginning  of  the  current  time  step.  KE  is  therefore  constant  for 
the  current  time  step. 


(2)  Is.  -  T*ySTiS 

It  is  generally  tedious  to  linearize  a  fourth  order  term  such  as  the  radiation 
heat  transfer  term  shown  above.  In  addition,  the  magnitude  of  the  radiation 
term  was  found  to  be  insignificant  when  compared  to  the  other  heat  transfer 
terms.  It  has  thus  been  neglected  for  the  course  of  this  work. 


(3)  fse  h(Te  —  Ta)NTdS 

Js  h(Te  —  Ta)NTdS  =  Js  hNTdSTe  -  ^  hNTTsdS 


=  [  hNTdSTe  -  [  /*NTN£dSTs  (3J1) 
J  Se  J  Se 


=  H eTe  -  KBTs 
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where  He  =  JSe  hNxdS,  and  K#  =  JSe  /iN^NydS. 


(4)  fSj  hlub{Td  -  Tw)NTdS 


Similar  to  the  procedure  used  in  (3),  this  term  can  be  expressed  as, 


* 


f  hlub{Td-Tw)NTdS  =  HdTd-KwTt 

J  St 


(3.12)  ► 


where  =  JSj  hiubNTdS ,  and  =  JSj  hiubNT~N^dS . 


(5)  Ssj  qfNrdS 

The  heat  generated  due  to  friction,  <jy,  is  calculated  as, 

9/  =  ffWr  | 

where  ff  is  the  frictional  stress,  and  from  Eq.  (2.25), 

ft  =  — mfcNT-^- 

|«r| 

Now,  we  have, 


f  qfNTdS  =  f  -mA:NrNjdSvr 

JSj  JSj 


=  [  -mfcNTNjdSvT  +  [  mkTSlTrfTxdSVd 
J  Sj  JSj 


=  — K  fVT  +  ~Bj)TVd 


(3.13) 


where  Ky  =  Js  mfcNrNjdS,  and  B dt  —  Ssj  mfcNj'NjTrdS.  The  explanation 
of  some  of  the  terms  described  above  may  be  obtained  from  section  2.3  describing 


* 
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the  interface  frictional  stresses.  Substituting  Eqs.(3.10)  through  (3.13)  into  Eq. 
(3.9),  and  expanding,  we  have, 

*  Qt  =  Kev  -  K fwT  -  KbTs  -  KttTtt  +  H dTd  +  H eTe  +  BDTVd  (3.14) 

Also,  from  Eq.  (3.8),  we  have, 

~  =  — Ct_1KcT  +  C  t_1Qt  (3-15) 

at 

Substituting  Eq.  (3.14)  into  Eq.  (3.15)  and  rearranging  the  terms,  the  full  size 
state  space  model  for  thermal  analysis  is  obtained  as, 

— —  =  A  dtV  +  AyT  +  BdtV<i  +  Wy  (3.16) 

at 

where  consists  of  the  matrices  Kjj,  K f,  and  Cy;  V  represents  the  velocity 

field  explained  in  section  2.4;  A t  is  built  from  Ks,  Km,  Kc  and  Cy;  and  Wy 
consists  of  H dTd  and  H eTe.  Similar  to  the  deformation  state  space  model,  all  the 
matrices  in  Eq.  (3.16)  are  built  at  time  step  t,  and  are  assumed  to  be  constant 
for  the  time  interval  At. 

3.3  Condensation  of  the  thermal  state  space  model 

In  nonisothermal  forging  analysis,  the  analytical  (state  space)  model  repre¬ 
sents  a  coupled  deformation-thermal  system  which  has  both  nodal  velocities  and 
nodal  temperatures  as  the  nodal  degrees  of  freedom.  The  total  number  of  finite 
'  element  equations  to  solve  is  (comparatively)  much  higher  than  in  isothermal 

analysis.  It  is  thus  necessary  to  condense  the  thermal  state  space  system  to  re- 
►  duce  the  number  of  degrees  of  freedom  involved.  Eq.  (3.16)  has  two  distinct 
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parts:  Aj^jV,  related  to  deformation,  and  AjT,  related  to  temperature  effects. 
Therefore,  system  condensation  is  performed  in  two  stages,  one  related  to  velocity 
field  condensation  and  the  other  dealing  with  temperature  field  condensation. 


* 


(1)  Condensation  of  the  velocity  field: 

Similar  to  the  condensation  of  the  deformation  model,  the  velocity  field  is 
partitioned  into  five  zones,  and  A dt  is  expanded  based  on  these  zones  as, 

A  dtV  =  AdtiVt  +  Adt2^n  +  AdtsVi  +  Adt-jVb  +  AdtsVc  (3-17) 

Now,  from  Eq.  (2.16),  we  have, 

Ksv  +  KtAv  =  F 


Because  Av  is  small  compared  to  v,  and  because  A/xrV  is  not  a  significant  part 
of  the  heat  flux  in  a  hydraulic  press  forging  process,  the  term  K* Av  is  left  out  to 
simplify  the  thermal  condensation  process.  Now,  we  can  solve  for  V £  from  the 
above  two  equations  as, 

VB  =  K-!4  [-K,4iVt  -  Ki43V7  -  KS42Tj,Vd]  (3.18) 


In  addition, 


Vjv  =  TyVd 


Substituting  V n  and  V £  into  Eq.  (3.17)  and  setting  Vc  =  0  (symmetric  bound¬ 
ary  condition),  we  have, 


A  otV  =  ApTl^T  +  ■^■DT2^I  +  B '  DT^d 


where, 


Adti  =  Apti  -  Az>T4Kj4iKa4i 


4 
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► 


A-DT2  —  A-DT3  ~  A£)T4K344Kj43 
B ' DT  =  (AdT2  —  Aot4Kj414Kj42)  Ty 
Substituting  these  expressions  into  Eq.  (3.16)  and  setting 

A  DT  =  [a -DTI  A-DT2]  B£T  =  DT  +  B£>T, 
the  thermal  state  space  model  after  velocity  condensation  may  be  represented 


as: 

—  =  ADr[y,Tl +ATT  +  BOTKi  +  Wr  (3'19) 

dt  LV/J 


Here,  the  retained  velocity  field  consists  of  nodal  velocities  of  the  die-contacting 
boundary  nodes  and  nodal  velocities  of  the  nodes  constituting  the  element  of  in¬ 
terest. 


(2)  Condensation  of  the  temperature  field: 

The  second  stage  of  condensation  is  done  with  respect  to  the  temperature 
field  so  as  to  retain  the  ‘critical’  nodal  temperatures.  To  be  consistent  with  the 
deformation  condensation,  the  ‘critical’  nodal  temperatures  are  identified  as  tem¬ 
peratures  of  the  nodes  at  the  die-workpiece  boundary,  and  temperatures  of  the 
nodes  constituting  the  element  of  interest.  In  accordance,  the  temperature  field 
is  divided  into  two  zones.  The  first  one  consists  of  the  temperatures  of  interest 
Tj,  which  in  turn  includes  the  nodal  temperatures  at  the  die-billet  interface  and 
the  nodal  temperatures  of  the  element(s)  of  interest.  The  rest  of  the  nodal  tem¬ 
peratures  are  put  into  the  second  group  represented  as  T m-  Eq.  (3.19)  is  now 


partitioned  and  rearranged  as  follows: 

d  Tj  _  Adtu  Adtu]  \vt]  +  Atu  Atu 
Jt[TM\~  [AdT2\  A-DT22  IT/]  Ut21  Aj22 . 


vd  + 


WTi‘ 

WT2. 

(3.20) 
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Assume  T m  is  the  unknown  temperature  field  at  (t  + At)  and  T^f  constitutes  the 
known  temperatures  at  time  f.  Using  the  forward  difference  scheme,  we  have, 

_  T m  ~  Tju  » 

dt  At 

Substituting  the  above  expression  into  the  2nd  row  of  Eq.  (3.20)  and  solving  for 
T m,  we  get, 

T M  =  At[I  —  At  Ay22]  1  [A-DT2iVt  +  A_D7'22V/  -i-  A^jT/  +  B DT2^d.  +  Wj2] 

+  [I  -  AtAT22]_1TM 

(3.21) 

where  I  is  the  unit  (or  identity)  matrix.  By  substituting  Eq.  (3.21)  into  the  1st 
row  of  Eq.  (3.20)  and  rearranging,  the  following  equation  is  obtained: 

yT  ]  +  KtpTi  +  B  TPVd  +  W  TP  (3.22) 

where, 

Kpr  =  [Kdti  K.DT2] 

Kuti  =  Adtu  +  ATl2Af[I  -  At Aj,22]  1Adt2  1 
Kdt2  —  A  qt\2  +  Aj12At[I  —  AtAj’22]  1A.£)T22 
K  tp  =  Ayu  +  Ay12At[I  —  AtAy22]  1Aj,2i 
Btp  =  BjDti  -l-  Ati2  At  [I  —  A£At22]  1Bdt2 
Wyp  =  Wn  +  Ayi2  At  [I  —  AtAy22]  1  Wj2  +  Aji2  [I  —  At  Aj22]  1Tj V/ 

Eq.  (3.22)  represents  the  condensed  state  space  model  for  thermal  analysis.  Here  ^ 
again,  all  the  matrices  are  built  at  time  t  and  are  considered  to  be  invariant  for 
the  time  interval  At.  -4 
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CHAPTER  4 


Coupled  Thermo-mechanical  State  Space  Model 


4.1  Construction  of  the  control  output  matrix 

As  mentioned  earlier,  for  this  work  the  effective  strain  rate  of  the  element  of 
interest  has  been  identified  as  an  important  process  variable  to  be  monitored  and 
controlled.  At  the  same  time,  since  nodal  velocities  (explained  earlier)  have  been 
chosen  as  the  state  variables,  it  is  necessary  to  formulate  an  output  equation 
to  represent  the  strain  rate  in  terms  of  the  nodal  velocities.  The  relationship 
between  the  strain  rate  and  nodal  velocities  is  a  nonlinear  function  over  the  entire 
domain  of  the  process.  Therefore,  at  every  simulation  time  step  a  linearization 
is  performed  to  formulate  the  output  equation. 

Let  v  represent  the  nodal  velocities  of  the  nodes  constituting  the  element  of 
interest.  Also  assume  the  current  state  space  model  is  built  with  the  velocity 
field  v*  and  is  valid  until  time  t.  Then,  in  the  next  time  step  (t  to  t  +  At),  the 
velocity  field  may  be  described  as, 

v  =  v*  +  Av  (4-1) 


From  Eq.  (2.18)  and  Eq.  (4.1), 

I2  =  vTPv 

=  (v*  +  Av)TP(v*  +  Av)  (4-2) 

=  v?Pv*  +  AvtPv*  +  v^PAv  +  AvtPAv 
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Because  the  last  term  in  the  above  equation  is  a  higher  order  term,  it  may  be 
neglected.  The  second  term  can  be  expanded  as, 


AvtPv»  =  AviPij j  v4 


=  £  Y.v*ipa  Avi 


*  \  3 

=  vJPTAv 


(4.3) 


From  Chapter  2,  P  is  equal  to  BjDB3,  and  D  is  a  diagonal  matrix  implying 
Pr  =  P.  Now,  from  Eq.  (4.3), 


AvTPv*  =vJPTAv 
-  vJPAv 

By  substituting  these  relations  into  Eq.  (4.2), 
t2  =  vfPv,  +  2v^PAv 

=  vJPv*  +  2v^PAv  +  vJ’Pv*  -  v^Pv* 
=  2v^Pv*  +  2v^PAv  -  v^Pv* 

=  2v;fP  (v*  +  Av)  -  vJPv* 

=  2v^Pv  —  vJPv, 


(4.4) 


(4.5) 


Now,  by  defining, 

2vJP  =  ci  vfPv,  =  ei 


a  hnearized  relationship  between  effective  strain  rate  and  nodal  velocities  is  ob¬ 
tained,  given  by, 

(f)2  =  civ-ei  (4.6) 


* 


v 


A 
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4.2  Coupled  state  space  system 


% 


A 


* 


► 


In.  the  state  space  model  for  deformation  analysis,  the  state  vector  is  given  by 
[  VT  ]T-  ft  the  strain-rate  squared  of  the  element  of  interest  (element  I)  is 
defined  as  tj,  the  following  output  equation  for  deformation  analysis  is  obtained: 


a2 

£I 


Ci 


VT 
L  V/ 


4-  Ei 


(4.7) 


Strain-rate  squared  is  used  as  the  output  variable  instead  of  the  strain-rate  be¬ 
cause  it  is  difficult  to  obtain  a  linearized  relationship  between  e  and  v.  Also, 
since  e  is  always  positive,  E^  gives  almost  the  same  information  as  E.  The  ci 
and  ei  in  Eq.  (4.6)  have  been  placed  in  the  corresponding  rows  and  columns 
of  Ci  and  Ei,  which  are  in  turn  related  to  the  nodal  velocities  associated  with 
the  element  of  interest.  In  the  thermal  state  space  model,  the  state  vector  T/ 
contains  the  nodal  temperatures  of  the  die- contacting  boundary  nodes,  and  the 
nodal  temperatures  of  the  element(s)  of  interest.  Theoretically  (and  from  an 
implementation  point  of  view),  it  is  difficult  to  include  all  of  these  temperatures 
in  the  optimal  design.  One  effective  way  to  circumvent  this  problem  is  to  select 
one  critical  nodal  temperature  from  the  set  Tj  to  be  the  output  variable.  For 
example,  if  the  objective  is  to  keep  all  the  nodal  temperatures  above  a  certain 
value,  then  the  lowest  nodal  temperature  could  be  selected  as  the  critical  nodal 
temperature.  With  the  critical  nodal  temperature  controlled,  the  rest  of  nodal 
temperatures  would  also  be  in  the  desired  temperature  range.  If  the  critical 
temperature  is  Tcj  the  output  equation  for  the  thermal  analysis  becomes, 


Tc  -  C2T / 


(4.8) 
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where,  C2  is  a  matrix  with  the  diagonal  element(s)  equal  to  1,  corresponding  to 
the  critical  nodal  temperature,  while  the  rest  of  the  diagonal  elements  are  zeros. 


In  this  study,  the  state  vector  is  defined  as  x,  where, 

XT  =  [  vj  Tj  ]  (4.9) 

and  V 5  and  T j  are  the  state  variables  corresponding  to  deformation  and  thermal 
analysis,  respectively.  V5  includes  the  nodal  velocities  of  the  die-contacting 
boundary  nodes  and  the  nodal  velocities  of  the  nodes  constituting  the  element 
of  interest.  T j  contains  the  nodal  temperatures  of  the  boundary  nodes  and  the 
nodal  temperatures  of  the  nodes  forming  the  element  of  interest. 

The  output  vector  for  the  coupled  thermo-mechanical  system  is  defined  as  y , 
where, 

yT  =  [  i)  T%  ]  (4.10) 

Here,  tj  is  the  strain  rate  of  the  element  of  interest  (explained  earlier),  and  Tc 
is  the  critical  nodal  temperature  as  discussed  above. 


With  the  above  definitions,  a  coupled  thermo-mechanical  state  space  model 
for  nonisothermal  metal  forming  processes  can  be  obtained  by  grouping  together 
the  deformation  state  space  model  of  Eq.  (2.38)  and  the  thermal  state  space 
model  of  Eq.  (3.22).  Then, 


x  =  Ax  +  BVrf  +  W 
y  -  Cx  4-  E 


(4.11) 


* 


y 


K 
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where  A  is  the  plant  matrix  and  is  given  by, 


% 

A  =  [  Ad  _0 

LK  dt  K  TP. 

B  is  the  input  matrix  given  by, 

< 

Bt  = 

[  B  l  B ttp  ] 

W  is  a  constant  vector  given  by, 

WT  =  [ 

f - 1 

E-<& 

S-hQ 

£ 

C  is  the  output  matrix,  where, 

n  _ 

rc,  o  1 

L  o  c2J 

and  E  is  a  constant  vector  in  (t,  t  +  At)  given  by, 

Et  =  [  Ef  0T  ] 

It  should  be  noted  that  W  contains  certain  terms  coming  from  the  condensa¬ 
tion  procedure.  In  keeping  these  terms,  the  error  between  the  condensed  model 
and  the  full  size  model  is  reduced  considerably,  especially  for  large  scale  problem- 
s.  Comparing  with  standard  state  space  models  generally  used  in  control  system 
theory,  Eq.  (4.6)  may  be  treated  as  a  system  with  constant  perturbation  W.  The 
model  is  built  and  updated  at  each  time  step  of  the  simulation.  An  unique  feature 
of  this  system  is  that  the  geometry  of  the  workpiece  and  the  boundary  condition 
keeps  changing  with  progress  in  deformation.  Because  of  this,  the  dimension  of 
the  plant  matrix  and  other  matrices  representing  the  system  also  changes  with 
increasing  stroke. 
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CHAPTER  5 


Process  Parameter  Design  Using  Optimal  Control  Theory 


5.1  Design  strategy 

To  establish  an  effective  design  strategy,  it  is  important  to  understand  the 
essential  characteristics  of  the  system  and  the  relationship  between  the  vari¬ 
ous  interacting  parameters  and  field  variables.  Chapters  2,  3  and  4  described  the 
different  field  variables  and  design  parameters  being  considered,  and  the  relation¬ 
ship  between  them.  This  section  describes  the  setting  up  of  the  optimal  control 
problem  for  forging  processes  and  the  general  strategy  adopted  for  solving  it. 

In  this  work,  the  design  goal  is  to  monitor  and  control  two  important  pro¬ 
cess  variables  associated  with  the  deformation  process,  namely  effective  strain 
rate  and  nodal  temperature.  Effective  strain-rate  is  directly  dependent  on  nodal 
velocities;  it  is  therefore  an  instantaneous  quantity  and  may  be  directly  influ¬ 
enced  by  changes  in  the  die  velocity.  However,  there  are  some  system  dependent 
constraints  on  how  effectively  the  strain  rate  may  be  controlled  while  simulating 
the  forging  process.  The  metal  working  process,  though  nonlinear  in  nature,  is 
treated  in  a  piecewise  continuous  manner  by  linearizing  it  at  every  simulation 
time  step.  The  deformation  analysis  therefore  does  not  allow  a  large  change  in 
effective  strain  in  one  time  step  of  the  simulation.  This  constraint  is  represented 
by  the  maximum  strain  increment  per  step  in  ALPID.  In  addition,  there  are 
other  constraints  such  as  the  forging  machine’s  acceleration  capability.  But  it  is  * 
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an  accepted  fact  in  the  metal  forming  field  that  strain  rate  in  the  billet  may  be 
controlled  by  suitably  adjusting  the  die  velocity.  Therefore,  a  design  methodolo¬ 
gy  based  on  control  theory  is  developed  to  design  optimal  ram  velocity  schedules 
that  meet  the  design  requirements  on  the  strain  rate.  Generally,  the  design  ob¬ 
jective  is  to  maintain  the  effective  strain-rate  at  some  desired  value,  or  within  a 
range  of  desired  values.  Considering  that  the  strain-rate  at  time  t  may  be  far 
away  from  the  desired  value,  to  achieve  the  above  objective  a  tracking  problem 
needs  to  be  formulated,  wherein  the  desired  value  may  be  reached  after  a  finite 
number  of  steps,  and  from  then  on  maintained  at  this  value.  At  each  time  step, 
an  optimal  die  velocity  is  generated  based  on  the  desired  strain-rate  value  at  that 
step,  which  is  defined  as  zn ,  where  4i’  stands  for  the  current  step  number.  Z\i  is 
decided  by  the  final  desired  strain-rate  value  (Id),  the  effective  strain-rate  value 
from  the  previous  step  (e{— 1),  and  the  specified  constraints.  This  may  be  stated 
mathematically  as, 


_  f  ki- 1  ±  mx£i- 1,  \zu  -  %d\  >  0.03td 

Zli"  lid,  |*ii-td|<0.03ld 


(5.1) 


where  mx  is  a  suitable  percentage  value  (usually  less  than  5%)  based  on  the 
machine’s  acceleration  capability  and  the  strain  increment  limit  [58]. 


Temperature  design  is  a  much  more  difficult  and  complicated  issue  than  veloc¬ 
ity  design.  From  chapter  3  it  is  observed  that  nodal  temperatures  are  influenced 
by  a  number  of  heat  generation  and  heat  transfer  terms.  The  first  term  in  Eq. 
(3.9)  (deformation  heat  generation)  is  dependent  on  the  type  of  forging  process 
used.  This  term  plays  a  more  important  role  in  high  speed  forging  processes  like 
mechanical  press  forging  (about  20  in/sec  ram  velocity)  or  hammer  (about  80 
in/sec  ram  velocity)  forging.  For  these  processes  changing  the  die  velocity  in  the 
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middle  of  the  process  is  virtually  impossible  because  of  the  high  inertia  of  the  die 
and  the  very  short  processing  time  involved.  Therefore,  in  this  study  we  confine 
ourselves  to  conventional  forging  processes  utilizing  hydraulic  presses  having  ram 
velocities  below  5  in/sec.  Earlier,  Boer  et.  al.  [59]  designed  a  varying  die  veloc¬ 
ity  profile  for  processing  NIM80A  material  and  implemented  it  on  a  hydraulic 
machine  with  a  microcomputer  as  the  controller.  Since  hydraulic  presses  oper¬ 
ate  at  relatively  low  die  velocities,  the  deformation  heat  generation  is  not  very 
significant.  Generally,  while  performing  forging  using  hydraulic- type  machines, 
the  temperature  in  the  middle  of  the  billet  increases  by  30  to  40°  E,  while  on 
the  boundary  surface  this  effect  is  totally  counteracted  by  the  die  chilling  effect, 
which  is  explained  below. 

The  last  four  terms  in  Eq.  (3.9)  mainly  affect  the  boundary  surface  tempera¬ 
ture.  Among  these  terms,  the  radiation  heat  term  is  neglected  because  it  has  an 
insignificant  magnitude  when  compared  to  the  other  terms.  Also,  if  the  process 
does  not  last  very  long,  the  convection  heat  loss  is  again  not  very  significant.  At 
the  boundaries,  the  friction  heat  has  some  influence  on  the  temperature,  but  is 
still  not  the  dominant  term.  The  dominant  factor  is  the  heat  conducted  between 
the  billet  and  die  through  the  interface.  Because  the  die  temperature  is  usually 
much  lower  than  that  of  the  billet  (up  to  1000°  F  at  times)  the  heat  transfer  due 
to  conduction  between  these  two  bodies  is  very  large,  and  a  severe  temperature 
gradient  exists  at  the  die-workpiece  boundary.  This  phenomenon  is  called  the 
die  chilling  effect.  Compared  to  the  40°  F  increase  in  temperature  in  the  middle 
of  the  billet  due  to  deformation  heat  generation,  the  boundary  temperature  has 
a  300  to  400° F  reduction  in  temperature  due  to  the  die  chilling  effect.  Due  to 
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the  large  temperature  gradient  between  the  center  of  the  billet  and  the  contact 
boundary,  severe  thermal  stresses  are  induced  which  can  radically  increase  die 
wear.  The  die  chilling  effect  also  affects  the  microstructure  and  mechanical  prop¬ 
erties  of  the  workpiece  because  the  resulting  temperature  gradient  could  result 
in  a  non-uniform  distribution  of  service  properties  in  the  billet.  Therefore,  in 
this  study,  one  of  the  primary  goals  is  to  determine  a  method  to  alleviate  the  die 
chilling  effect.  Since  this  phenomenon  occurs  mainly  at  the  boundaries,  only  the 
treatment  of  boundary  nodal  temperatures  was  taken  up  during  this  work. 

The  difficulties  involved  in  boundary  nodal  temperature  control  are  appar¬ 
ent  from  its  following  two  characteristics.  First,  temperature  is  an  accumulated 
value  and  not  an  instantaneous  quantity  like  strain  rate.  Second,  it  is  not  di¬ 
rectly  influenced  by  the  die  velocity.  So,  changes  made  in  die  velocity  at  some 
instant  in  time  do  not  immediately  affect  the  temperature.  In  this  work,  the  ob¬ 
jective  selected  for  temperature  control  is  to  reduce  the  temperature  range  in  the 
workpiece  material.  This  is  desirable  because  it  would  result  in  a  more  uniform 
distribution  of  mechanical  properties  in  the  final  product.  Since  temperature  is 
not  easily  controllable,  for  the  scope  of  this  work  the  goal  is  confined  to  raising 
the  nodal  temperatures  above  a  specified  value  rather  than  trying  to  maintain 
them  at  the  given  value.  To  achieve  the  above  objectives,  the  design  scheme  is 
again  formulated  as  a  tracking  problem,  and  the  lowest  boundary  nodal  temper¬ 
ature  is  selected  as  the  design  parameter  to  achieve  this  goal.  The  logic  behind 
this  is  if  the  lowest  boundary  temperature  is  maintained  above  a  certain  value,  all 
the  other  boundary  nodal  temperatures  would  also  be  above  the  specified  value. 
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Similar  to  strain-rate  control,  we  define  the  desired  nodal  temperature  at  step 
i  as  Z2i,  which  can  be  generated  in  two  ways.  One  way  is  to  use  the  following 
formula: 

Z2i  =  Ti_i  -f  gdi-i  (1  —  n)  At  (5.2) 

where  T{-\  is  the  nodal  temperature  at  step  (i-1 ),  and  gd{-i  is  the  temperature  ^ 

gradient  at  step  (i-1).  n  is  a  positive  percentage  value  usually  less  than  3%  and 

is  again  determined  based  on  the  machine  capabilities  and  the  strain  increment 

limit.  The  objective  of  Eq.  (5.2)  is  to  reduce  the  temperature  range  by  reducing 

the  temperature  gradient  at  each  simulation  step.  Though  it  cannot  give  an 

accurate  control  of  temperature,  it  is  quite  suitable  for  nonisothermal  design 

involving  die  velocity  as  the  lone  design  variable.  The  other  method  of  defining 

the  desired  nodal  temperature  is  by  inputting  a  known  polynomial  approximation 

of  the  required  temperature  into  the  design  scheme,  as  discussed  later  in  the 

report. 

In  the  nonisothermal  design  process,  the  output  vector  contains  at  least  two 
output  variables:  the  effective  strain-rate  and  the  nodal  temperature.  Without 
the  initial  die  temperature  parameter  (explained  later),  the  only  control  input 
is  the  die  velocity.  From  control  system  theory,  this  is  a  single  input  multiple 
output  (SIMO)  system.  The  forging  system  also  contains  a  constant  vector  W  in 
the  system  equation,  which  may  be  treated  as  a  constant  perturbation  term  [60]. 

For  problems  of  this  nature,  several  kinds  of  optimal  control  schemes  have  been 
developed.  For  example,  the  robust  servomechanism  design  in  reference  61  uses  1 

two  servocompensators  so  that  the  resulting  system  is  stable,  and  the  steady-state 
error  is  zero  for  all  kinds  of  disturbances.  That  means  the  disturbance  could  occur  ^ 
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in  the  system  matrices,  input /output  vectors,  or  the  feedback  channels.  The 
Proportional  plus  Integral  state  feedback  design  in  reference  60  uses  an  integral 
feedback  to  eliminate  the  steady-state  error  caused  by  a  known  or  unknown 
constant  disturbance.  However,  the  above  two  design  schemes  have  the  basic 
requirement  that  the  number  of  control  inputs  should  not  be  less  than  the  number 
of  system  outputs.  Therefore,  this  scheme  cannot  be  utilized  in  the  current  SIMO 
problem. 

In  this  work,  the  linear  quadratic  regulator  (LQR)  theory  is  used  as  a  design 
tool  for  determining  the  optimal  (forging)  process  parameters.  LQR  theory  is 
one  of  the  most  popular  design  methods  in  optimal  control.  It  is  a  proportional 
feedback  design  and  can  be  used  for  SIMO  and  SISO  systems.  In  this  study,  a 
terminal  condition  is  also  introduced  into  the  standard  LQR  quadratic  perfor¬ 
mance  index  to  improve  the  error  characteristics  of  the  system.  From  Eq.  (4.11), 
for  step  i  which  lasts  from  time  t  to  (t  +  At),  the  desired  value  vector  may  be 
defined  as, 

z i  =  [zii  Z2i\ 

and  the  error  vector  may  be  expressed  as, 


ei(t)  =  Zi-yi(t)  te{t,t  +  At) 

Now,  the  following  performance  index  may  be  used  in  the  optimal  design  scheme: 

J  =  1  e?  (At)  Me;  (At)  +  \  jf"  [ef  (f)  Qe,  (t)  +  Vj  (f)  RUf  (t)]  it  (5.3) 

where  M,  Q,  and  R  are  the  weighting  matrices  for  the  terminal  condition,  error, 
and  input  term,  respectively,  and  U{  is  the  input  vector.  For  the  isothermal  case, 
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or  in  cases  where  only  die  velocity  design  is  considered,  U{(t)  =  [V*-(t)].  If  initial 
die  temperature  design  is  also  performed,  U J(t)  =  \Ydi{^)  AT^f)],  where  A T^i 
is  the  initial  die  temperature  adjustment  parameter.  Minimization  of  Eq.  (5.3)  y 

gives  [60]: 

U i  (f )  =  — R-1Bt  [Kt  (t)  x  (i)  -  g i  (<)]  (5.4) 

where  Ki(t)  and  g{(t)  are  the  solutions  to  the  differential  Riccati  equations  given 

by, 

Ki  (t)  =  -Ki  (0  A  -  ATKt  (t)  +  Ki  ( t )  BR-1BTKi  (t)  -  CTQC 

g i  (t)  =  [Ki  (t)  BR-1Bt  -  At]  gi  ( t )  -  CtQz i  +  Ki  (0  Wr  (5.5) 

In  Eq.  (5.5),  Ki  and  gi  are  time  dependent,  and  are  not  constant  because 
of  the  finite  time  integral  in  the  performance  index.  However,  if  the  output 
is  more  than  the  control  input,  there  may  be  no  solution  for  the  infinite  time 
optimal  tracking  problem  [60].  Eq.  (5.4)  gives  a  solution  series  in  the  time 
interval  (f,f  +  At).  To  make  the  process  more  efficient  and  for  reducing  the 
computational  time,  the  average  of  this  solution  series  is  taken  as  the  optimal 
design,  i.e., 

XJ°pt  =  -1 -  tj  €  (i,  t  +  At )  (5.6) 
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where  n  is  the  number  of  smaller  time  intervals  into  which  step  i  is  divided  for 
solving  Eq.  (5.5).  For  most  cases,  it  has  been  observed  that  Eq.  (5.6)  is  a 
reasonably  good  and  acceptable  assumption.  y 

To  give  an  overview  of  the  methodology,  the  state  space  model  is  developed 

4 

based  on  the  finite  element  simulation,  and  then  the  LQR  design  scheme  is  used 
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to  design  the  optimal  process  parameters.  This  is  an  off-line  design  scheme,  and 
the  purpose  is  to  use  control  theory  as  a  design  tool  rather  than  to  design  a 
practical  control  system.  The  emphasis  is  on  obtaining  die  velocity  profiles  and 
the  initial  die  temperature,  rather  than  in  stabilizing  the  control  loops  as  done  in 
classical  control  system  design.  The  entire  design  scheme  developed  is  illustrated 
in  Fig.  6. 

During  nonisothermal  forging  process  design,  the  strain-rate  control  and  the 
temperature  control  may  have  a  conflicting  requirement  on  the  die  velocity.  For 
example,  the  strain  rate  control  may  need  the  die  velocity  to  be  reduced;  at  the 
same  time  the  temperature  control  may  actually  require  the  die  velocity  to  be 
increased.  In  such  a  situation  the  strain-rate  control  is  given  higher  priority  be¬ 
cause  the  strain-rate  is  more  sensitive  to  changes  in  the  die  velocity.  Therefore, 
in  this  design  scheme,  if  the  above  mentioned  conflict  occurs,  the  die  velocity 
is  decided  based  on  the  strain-rate.  Actually,  temperature  is  not  an  easily  con¬ 
trollable  quantity  and  the  main  contribution  of  the  optimal  die  velocity  design 
to  the  thermal  aspect  of  the  process  is  in  reducing  the  temperature  range  in 
the  deforming  billet  by  reducing  the  die-workpiece  contact  time.  This  may  not 
always  be  possible,  though,  and  a  more  effective  way  to  directly  influence  the 
interface  temperature  is  by  designing  an  appropriate  initial  die  temperature  for 
the  process. 


t 


4 


1 


Figure  6:  Design  scheme  flow  chart 
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5.2  Initial  die  temperature  design 


4 


Because  of  the  strong  die  chilling  effect  during  hydraulic  press  forging,  the 
most  efficient  way  to  influence  the  boundary  nodal  temperatures  and  reduce  the 
temperature  range  in  the  workpiece  is  by  designing  the  initial  die  temperature. 
Since  the  entire  design  methodology  is  based  on  the  LQR  design  scheme,  it  is 
preferable  to  continue  to  apply  the  LQR  method  for  initial  die  temperature  design 
also.  But,  LQR  is  a  step  by  step  design  scheme,  while  initial  die  temperature 
is  a  single  value  and  needs  to  be  specified  at  the  beginnning  of  the  process.  To 
solve  this  problem,  the  following  design  approach  is  proposed. 


In  this  scheme,  a  new  design  variable  called  die  temperature  adjustment 
(A Td)  is  introduced.  By  adding  AT^  to  the  present  die  temperature  ( Td )  the 
total  die  temperature  is  given  by  ( Td  +  AT*),  and  Eq.  (3.12)  may  now  be  written 
as, 

[  hlub  (! Td  +  A Td  -  Tw)  N TdS  =  +  HdATd  -  KWTW 

JSi 


By  using  this  equation  instead  of  Eq.  (3.12)  in  developing  the  state  space  model, 
the  modified  state  equation  may  be  written  as, 


x  =  Ax  + 


bd 

0  1 

^  1 

.B  TP 

Hdj 

LatJ 

+  W  r 


(5.7) 


The  difference  between  Eq.  (4.11)  and  Eq.  (5.7)  is  that  in  the  latter,  the  input 
matrix  B  has  been  expanded  and  now  there  are  two  input  variables  instead  of 
one.  If  A Td  -  0,  then  Eq.  (5.7)  will  yield  the  same  results  as  Eq.  (4.11).  The 
term  H^AT^  actually  changes  the  temperature  gradient  at  each  step  by  a  certain 
amount.  mainly  contains  shape  functions  and  heat  transfer  coefficients  and 
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does  not  change  much  during  the  whole  process.  For  a  case  in  which  the  tem¬ 
perature  gradient  needs  to  be  reduced  by  the  same  amount  at  each  time  step  (by 
applying  the  A Td  correction),  nearly  the  same  A Td  value  may  be  expected  from 
the  LQR  design  scheme.  This  suggests  that  to  achieve  the  design  objective,  the 
actual  initial  die  temperature  should  be  changed  by  the  average  value  of  A Td  for 
all  the  time  steps,  i.e., 

EA  Td 

ATd0  =  - -  (5-8) 

n 

where  i  is  the  design  step  count  and  n  is  the  total  number  of  design  steps.  It 
has  been  found  that  with  the  (new)  initial  die  temperature  ( Td  +  A Td0)  designed 
in  this  manner,  the  boundary  nodal  temperature  may  be  altered  sufficiently  to 
meet  the  design  requirement. 

To  implement  this  scheme,  the  key  issue  during  initial  die  temperature  design 
is  how  to  generate  a  desired  nodal  temperature  profile  which  not  only  satisfies  the 
design  requirement,  but  also  gives  nearly  the  same  temperature  gradient  change 
at  each  design  step.  Through  parametric  studies  it  was  found  that  for  hydraulic 
press  forging,  the  temperature  curve  of  a  die- contacting  node  is  nearly  a  straight 
line.  Therefore,  we  may  with  reasonable  accuracy  use  a  3rd  order  polynomial 
to  approximate  the  temperature  path  of  a  particular  node  with  respect  to  time. 
Then,  by  changing  the  coefficient  of  the  first  order  term  in  the  polynomial,  we 
can  generate  the  desired  temperature  profile.  The  problem  now  becomes  one 
of  tracking  the  desired  temperature  profile.  The  above  design  scheme  works 
successfully  if  the  die  temperature  is  not  very  close  to  the  billet  temperature, 
and  if  the  die  velocity  does  not  vary  by  a  significant  amount.  These  conditions 
hold  good  in  most  hot  die  forging  cases.  The  detailed  procedure  for  designing 
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the  initial  die  temperature  is  explained  in  the  numerical  examples  section  of  this 
report. 

5.3  LQR  design  scheme 

In  Eq.  (4.11)  E  is  a  constant  vector  and  can  be  grouped  into  the  left  side 
of  the  output  equation,  with  the  vector  of  desired  values.  Now,  the  state  space 
model  representing  the  system  may  be  written  as, 

x  =  Ax  +  Bu  +  W  (5-9) 

y  =  Cx  (5.10) 

where  the  state  vector  x  has  a  dimension  n,  the  input  vector  u  has  a  dimension 
m  and  the  output  vector  y  has  a  dimension  p.  tl,  m,  and  p  refer  to  the  number 
of  states,  number  of  inputs,  and  number  of  outputs  in  the  system,  respectively. 
Then,  A  (system  matrix)  is  a  n  x  n  matrix,  B  (input  matrix)  isanxm  matrix, 
and  C  (output  matrix)  is  a  p  X  n  matrix.  W  is  a  known  vector  resulting  from 
the  finite  element  condensation  procedure.  During  the  derivation,  to  make  the 
notations  simple,  the  integration  limits  are  set  as  to  =  <jfc  and  T  =  tk  +  At, 
where  f0  and  T  are  the  lower  and  upper  integration  limits  in  the  performance 
index  equation.  If  the  desired  output  is  z,  and  the  error  vector  is  e  =  z  -  y,  the 
performance  index  is  given  by, 

J  =  leT  (T)  Me  (T)  +  \  F  [eT  (t)  Q  (t)  e  (<)  +  uT  (t)  Ru  (t)]  dt 

2  ^  J  to 

where  Q  is  a  n  xn  semi-positive  definite  matrix  for  weighting  the  error  vector, 
Risamxro  positive  definite  matrix  for  weighting  the  input  vector,  and  M  is 
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a  semi-positive  definite  matrix  used  for  weighting  the  terminal  condition.  The 
method  of  selection  of  these  weighting  matrices  is  explained  in  section  5.4  of  this 
chapter.  Some  important  steps  in  the  derivation  of  the  optimal  control  scheme 
are  described  below  [49].  The  error  vector,  e,  may  be  written  as, 


e  =  z  -  Cx  (5.11) 

By  following  the  standard  state  space  equation  solving  procedure,  and  by  sub¬ 
stituting  Eq.  (5.11)  into  the  performance  index,  the  Hamiltonian  function  may 
be  written  as, 

H  =  -[z  -  Cx]TQ  [z  —  Cx]  -f  ^utRu  -f  xTATA  +  uTBTA  +  WTA 

A  Z 

where  A  is  a  Lagrange  multiplier.  From  the  condition,  ^  =  0,  we  have, 

Ru  +  BtA  =  0 

u  =  -R-1BrA  (5.12) 

•  QTT 

From  A  =  —  -g^,  we  have, 

A  =  -CtQCx  +  CtQz  -  AtA  (5.13) 

By  substituting  Eq.  (5.12)  into  Eq.  (5.9), 

x  =  Ax  -  BR_1BtA  +  W  (5.14) 

Now,  by  defining  the  following  terms, 

s  =  br_1bt 

V  =  CrQC 


4 


I 

4 
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(5.15) 


P  =  CtQ 


we  can  group  Eq.  (5.13)  and  Eq.  (5.14)  as, 


A 


This  represents  a  2n  x  2ra  system,  where  n  is  the  number  of  states  in  the  state 
space  model. 


If  $  is  the  state  transition  matrix,  the  solution  to  Eq.  (5.15)  is  given  by, 


[*(T)1 

U(T)J 


W  (r)r(r)' 
.  p(r)z(T)  - 


(5.16) 


i.e. 


f*(T)l 

U(T)J 


$(T,t) 


x(t)' 

U(oJ 


+ 


gi  w 

-  &2  ( t ). 


(5.17) 


The  2 n  boundary  conditions  corresponding  to  the  current  system  are, 


t  =  to,  x(t0)  =  £ 

A(r’  =  3^f)[rT(r)Me(r) 


t  =  T; 


=  CT  (T)  MC  (T)  x  (T)  -  CT  (T)  Mz  (T) 


Now,  by  defining  Ct(T)MC(T)  =  G (T),  we  have, 


% 


A(T)  =  G(T)x(T)-Ct(T)Mz(T)  (5.18) 


Also,  say, 


$11  $12 
$2i  $12. 
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Then,  Eq.  (5.17)  becomes, 


fx(T)l 

#11  #12 

[x(T)l 

_L 

gi  (O' 

U(T)J 

#2i  $12 

UmJ 

1 

.  g2  (0 . 

Substituting  Eq.  (5.18)  into  the  above  equation,  we  have: 

A  (0  =  [#22  -  G#12]_1  ([G#n  -  #21]  X  ( t )  +  [Ggl  (t)  -  g2  (t)  -  CT  (T)  Mz]) 
Now,  by  assuming, 

K  (0  =  [#22  -  G#12]_1  [G#n  -  #21]  (5.19) 

g  (0  =  [#22  -  G#12]_1  [Ggl  ( t )  -  g2  (0  -  CT  (T)  Mz]  (5.20) 


we  have, 


A  =  K(t)x(t)-g  (t) 


(5.21) 


where  K (t)  is  a  n  x  n  time-varying  matrix,  and  g(£)  is  a  n- dimensional  vector. 

By  differentiating  Eq.  (5.21)  with  respect  to  t,  we  have, 

A  =  Kx  +  Kx  —  g  (5.22) 


Then,  substituting  Eq.  (5.21)  into  Eq.  (5.14),  we  get, 

x  =  [A  —  SK]  x  +  Sg  +  Wr 


(5.23) 


Now,  substituting  Eq.  (5.23)  into  Eq.  (5.22),  we  get, 

A  =  [K  +  KA  -  KSK]  x  +  KSg  +  KWr  -  g 


(5.24) 
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Substituting  Eq.  (5.21)  into  Eq.  (5.13),  A  may  also  be  written  as, 


A  =  [-V  -  ATKj  x  +  ATg  +  Pz  (5.25) 

\ 


A 


Comparing  Eq.  (5.24)  and  Eq.  (5.25),  we  have, 

K  =  -KA  -  AtK  +  KSK  -  V 


i.e., 

K  =  -KA  -  AtK  +  KBR"1BtK  -  CtQC 


(5.26) 


and, 

g=  [KBR“1BT-AT]g-CTQz  +  KWr'  (5.27) 


The  above  equations  are  the  differential  Riccati  equations  mentioned  in  section 
(5.1).  The  corresponding  boundary  condition  for  this  system  is: 


A(T)  =  K(T)x(T)-g(T) 

A  (T)  =  CT  (T)  MC  (T)  -  CTMz  (T) 

where 

K(T)  =  Ct(T)MC(T) 
g(T)  =  CT(T)Mz(T) 


(5.28) 


By  solving  Eq.  (5.26)  and  Eq.  (5.27)  together  with  the  boundary  conditions 
from  Eq.  (5.28),  the  optimal  design  solution  may  be  obtained.  While  computing 
the  solutions  of  Eq.  (5.26)  and  Eq.  (5.27),  the  Pade  approximation  [47]  method 
has  been  used  in  determining  the  state  transformation  matrix  $.  Eqs.  (5.19) 
and  (5.20)  are  then  used  to  calculate  K  and  g  [48]. 
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5.4  Selection  of  weighting  matrices  for  control 


The  selection  of  weighting  matrices  during  optimal  control  is  a  very  difficult 
issue,  and  no  general  method  exists  which  is  applicable  to  all  kinds  of  systems 
and  problems.  Generally,  weighting  matrices  are  determined  by  a  trial  and  error 
process  or  are  based  on  parametric  studies.  In  nonisothermal  forging  process 
design,  the  boundary  conditions  may  change  significantly  from  one  time  step 
to  the  next,  so  that  the  system  equations  could  be  totally  different  between 
steps.  Therefore,  using  constant  weighting  matrices  for  different  design  problems 
is  practically  impossible.  In  this  study,  the  following  matrices  are  used  during 
process  parameter  design  as  the  weighting  matrices  [65]: 

Q  =  diag 


1 

.  emax . 


M  =  100Q 


R  = 


On 


0 

0.01 


Here,  r  is  a  positive  value  varying  between  1  and  100,  and  is  automatically  decided 
by  the  computer  program  based  on  the  deviation  of  the  process  parameters  from 
their  respective  desired  values.  This  quantity  actually  forces  the  change  of  the 
design  parameters  within  fixed  constraints.  Usually,  an  isothermal  design  uses  a 
large  value  of  r  compared  to  a  nonisothermal  design  process.  For  plane  strain 
cases  in  general,  and  for  complex  problems  it  was  found  that  the  above  weighting 
matrix  scheme  was  resulting  in  numerical  difficulties  during  implementation.  For 
such  cases,  the  weighting  matrices  were  constrained  within  certain  finite  bounds 
determined  using  systematic  parametric  studies. 
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CHAPTER  6 


Numerical  Examples  -  Full  Size  State  Space  Model 

A  general  purpose  computer  program  COPP  (Control  of  Optimal  Processing 
Parameters)  has  been  developed  to  implement  the  design  techniques  described  in 
this  work.  This  program  is  built  on  the  rigid  viscoplastic  finite  element  program 
ALPID  (Analysis  of  Large  Plastic  Incremental  Deformation).  In  addition  to 
using  ALPID  subroutines  for  analysis  purposes,  a  design  and  control  module  has 
been  developed  and  built  into  the  program  for  designing  the  optimal  process 
parameters.  Along  with  the  ALPID  input  file,  COPP  requires  a  small  design 
input  data  file,  the  format  of  which  is  presented  in  Appendix  (B).  Appendix  B 
also  contains  a  detailed  user’s  manual  with  examples  for  running  this  program. 
In  general,  COPP  has  the  following  capabilities: 

1.  ALPID  simulation  for  both  isothermal  and  nonisothermal  processes. 

2.  Optimal  die  velocity  design  for  isothermal  process. 

3.  Optimal  die  velocity  design  for  nonisothermal  process. 

4.  Initial  die  temperature  design  for  nonisothermal  process. 
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This  chapter  presents  a  few  numerical  examples  to  support  and  substantiate 
the  theory  presented  in  earlier  sections  of  the  report.  The  case  of  an  axisymmetric 
engine  disk  forging  (Fig.  7)  is  considered  under  different  operating  conditions. 
The  material  used  for  the  simulation  is  Ti-6242  (Titanium  alloy),  which  is  a 
strain-rate  sensitive  material.  Due  to  its  symmetry,  only  a  quarter  section  of  the 
cylindrical  billet  is  modeled.  Discretization  of  the  billet  continuum  was  performed 
using  50  isoparametric  quadrilateral  elements  (Fig.  8),  and  the  die  was  modeled 
using  73  finite  elements.  A  constant  shear  friction  factor  of  0.3  was  used  at  the 
die-billet  interface.  A  billet  temperature  of  1735°F  and  a  die  temperature  of 
600° F  was  specified  at  the  beginning  of  the  process.  A  total  die  stroke  of  0.50 
inch  was  required  to  complete  the  forging  operation  as  shown  in  Fig.  7. 


6.1  Model  validation 

The  accuracy  of  the  state  space  model  in  representing  the  metal  forming 
system  has  been  tested  extensively  and  validated  using  ALPID  simulation  re¬ 
sults.  The  testing  was  done  for  both  isothermal  and  nonisothermal  processes 
using  varying  frictional  conditions,  different  die  velocities,  and  for  a  range  of  dif¬ 
ferent  materials.  The  response  of  the  state  space  system  was  found  to  be  in  close 
conformance  with  the  ALPID  simulation  results.  One  such  example  is  present¬ 
ed  here  to  establish  the  effectiveness  of  the  state  space  model  in  capturing  the  * 

deformation  characteristics  of  the  forging  process. 
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Figure  7:  Engine  disk  forging  -  FEM  simulation 
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Figure  8:  FEM  discretization  of  the  billet 
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Consider  the  forging  of  the  engine  disk  described  earlier.  Element  23  is 
selected  as  the  element  of  interest.  The  rationale  behind  this  is,  due  to  its  location 
in  the  workpiece,  element  23  is  always  under  deformation  and  is  also  likely  to 
go  through  a  number  of  modes  of  deformation.  As  a  result  this  element  is  very 
sensitive  to  changes  in  the  die  velocity  and  is  a  good  candidate  to  be  selected  as 
the  control  element.  As  shown  in  Fig.  9,  the  state  space  model  is  built  at  18% 
reduction  when  the  stroke  is  0.33  inch  and  the  die  velocity  is  1.0  in/sec.  In  the 
current  state  space  system,  17  state  variables  are  used  for  describing  the  nodal 
velocities.  Of  these,  9  are  the  tangential  nodal  velocities  of  the  die- contacting 
boundary  nodes  and  the  rest  are  nodal  velocities  of  element  23.  In  addition,  14 
states  are  used  to  describe  the  nodal  temperatures.  Among  these,  10  are  nodal 
temperatures  of  the  die-contacting  nodes,  and  the  other  four  are  temperatures 
of  the  nodes  constituting  element  23.  The  full  size  state  space  model  thus  has 
31  state  variables.  Now,  a  small  perturbation  of  0.0025  in/sec  is  given  to  the  die 
velocity  to  activate  the  response  of  the  state  space  model.  At  the  same  time,  this 
change  in  die  velocity  is  also  implemented  in  ALPID.  The  results  from  the  state 
space  system  and  those  from  ALPID  (which  serves  as  the  reference)  are  then 
tabulated  and  compared  (Table  1).  The  results  presented  in  Table  1  show  that 
the  temperature  field  and  the  nodal  velocity  field  (for  element  23)  generated  by 
the  state  space  system  closely  match  with  those  of  the  ALPID  simulation. 
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Figure  9:  Validation  of  the  state  space  model 
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Table  1.  Model  Validation:  State  space  model  vs.  ALPID 


t=0. 33292  +  0.0  (sec) 
Velocity  Field  (in/sec) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

23 

0.10078E+00 

0.10078E+00 

— 

35 

0.15708E+00 

0.15708E+00 

— 

47 

0.20641E+00 

0.20641E+00 

— 

59 

0.14646E+00 

0.14646E+00 

— 

71 

0.14127E-01 

0.14127E-01 

— 

125 

-0.11512E+00 

-0.11512E+00 

— 

127 

-0.26875E+00 

-0.26875E+00 

— 

129 

-0.29317E+00 

-0.29317E+00 

— 

131 

-0.33610E+00 

-0.33610E+00 

— 

53 

0.86712E+00 

0.86712E+00 

— 

54 

-0.30277E+00 

-0.30277E+00 

— 

55 

0.73772E+00 

0.73772E+00 

— 

56 

-0.52864E+00 

-0.52864E+00 

— 

65 

0.83599E+00 

0.83599E+00 

— 

66 

-0.14027E+00 

-0.14027E+00 

— 

67 

0.75102E+00 

0.75102E+00 

— 

68 

-0.21921E+00 

-0.21921E+00 

— 

N.O.F.:  Nodal  degree  of  freedom 
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Table  1-  (Continued) 


t=0. 33292  4-  0.0  (sec) 
Temperature  Field  ( °F ) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

6 

0.17418E+04 

0.17418E+04 

— 

12 

0.15857E+04 

0.15857E+04 

— 

18 

0.15509E+04 

0.15509E+04 

— 

24 

0.15602E+04 

0.15602E+04 

— 

27 

0.17465E+04 

0.17465E+04 

— 

28 

0.17491E+04 

0.17491E+04 

— 

30 

0.15454E+04 

0.15454E+04 

— 

33 

0.17421E+04 

0.17421E+04 

— 

34 

0.17433E+04 

0.17433E+04 

— 

36 

0.16389E+04 

0.16389E+04 

— 

63 

0.17079E+04 

0.17079E+04 

— 

64 

0.16738E+04 

0.16738E+04 

— 

65 

0.16513E+04 

0.16513E+04 

— 

66 

0.16337E+04 

0.16337E+04 

— 

.O.F.:  Nodal  degree  of  freedom 
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Table  1.  (Continued) 

t=0.33292  +  0.005  (sec) 
Velocity  Field  (in/sec) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

23 

0.10090E+00 

0.10562E+00 

4.4684 

35 

0.15719E+00 

0.16308E+00 

3.6095 

47 

0.20655E+00 

0.21551E+00 

4.1550 

59 

0.14656E+00 

0.15717E+00 

6.7457 

71 

0.14061E-01 

0.14151E-01 

0.6354 

125 

-0.11519E+00 

-0.12870E+00 

10.4985 

127 

-0.26891E+00 

-0.27671E+00 

2.8208 

129 

-0.29331E+00 

-0.29747E+00 

1.3996 

131 

-0.33622E+00 

-0.33604E+00 

0.0537 

53 

0.86752E+00 

0.88153E+00 

1.5902 

54 

-0.30291E+00 

-0.30300E+00 

0.0301 

55 

0.73804E+00 

0.75044E+00 

1.6529 

56 

-0.52889E+00 

-0.53073E+00 

0.3459 

65 

0.83638E+00 

0.84936E+00 

1.5283 

66 

0.14034E+00 

-0.13878E+00 

1.1298 

67 

0.75135E+00 

0.76265E+00 

1.4823 

68 

-0.21932E+00 

-0.21743E+00 

0.8694 

N.O.F.:  Nodal  degree  of  freedom 
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Table  1.  (Continued) 


t— 0.33292  +  0.005  (sec) 
Temperature  Field  (°  F) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

6 

0.17420E+04 

0.17420E+04 

0.0006 

12 

0.15840E+04 

0.15842E+04 

0.0103 

18 

0.15489E+04 

0.15488E+04 

0.0072 

24 

0.15584E+04 

0.15582E+04 

0.0096 

27 

0.17467E+04 

0.17466E+04 

0.0044 

28 

0.17494E+04 

0.17494E+04 

0.0039 

30 

0.15432E+04 

0.15430E+04 

0.0122 

33 

0.17423E+04 

0.17423E+04 

0.0032 

34 

0.17435E+04 

0.17434E+04 

0.0058 

36 

0.16370E+04 

0.16366E+04 

0.0235 

63 

0.17063E+04 

0.17060E+04 

0.0172 

64 

0.16713E+04 

0.16714E+04 

0.0023 

65 

0.16490E+04 

0.16490E+04 

0.0014 

66 

0.16312E+04 

0.16311E+04 

0.0047 

.O.F.:  Nodal  degree  of  freedom 
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Table  1.  (Continued) 


t=0. 33292  +  0.01  (sec) 
Velocity  Field  (in/sec) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

23 

0.10099E+00 

0.10842E+00 

6.8529 

35 

0.15728E+00 

0.16655E+00 

5.5713 

47 

0.20665E+00 

0.22100E+00 

6.4914 

59 

0.14664E+00 

0.16244E+00 

9.7278 

71 

0.14011E-01 

0.15191E-01 

7.7659 

125 

-0.11524E+00 

-0.13458E+00 

14.3705 

127 

-0.26902E+00 

-0.27961E+00 

3.7885 

129 

-0.29342E+00 

-0.29832E+00 

1.6414 

131 

-0.33634E+00 

-0.33486E+00 

0.4417 

53 

0.86784E+00 

0.88620E+00 

2.0716 

54 

-0.30302E+00 

-0.30269E+00 

0.1094 

55 

0.73830E+00 

0.75574E+00 

2.3080 

56 

-0.52910E+00 

-0.53030E+00 

0.2280 

65 

0.83671E+00 

0.85311E+00 

1.9227 

66 

-0.14040E+00 

-0.13830E+00 

1.5162 

67 

0.75161E+00 

0.76531E+00 

1.7895 

68 

-0.21940E+00 

-0.21693E+00 

1.1403 

N.O.F.:  Nodal  degree  of  freedom 


Table  1.  (Continued) 


t=0. 33292  +  0.01  (sec) 


Temperature  Field  (°  F) 


N.O.F. 

State  Space  Model 

ALPID 

Error  (%) 

6 

0.17423E+04 

0.17423E+04 

0.0010 

12 

0.15824E+04 

0.15827E+04 

0.0212 

18 

0.15469E+04 

0.15467E+04 

0.0138 

24 

0.15565E+04 

0.15562E+04 

0.0185 

27 

0.17470E+04 

0.17468E+04 

0.0088 

28 

0.17497E+04 

0.17496E+04 

0.0078 

30 

0.15409E+04 

0.15406E+04 

0.0232 

33 

0.17425E+04 

0.17424E+04 

0.0064 

34 

0.17437E+04 

0.17435E+04 

0.0116 

36 

0.16350E+04 

0.16343E+04 

0.0455 

63 

0.17047E+04 

0.17041E+04 

0.0373 

64 

0.16689E+04 

0.16690E+04 

0.0050 

65 

0.16466E+04 

0.16466E+04 

0.0031 

66 

0.16286E+04 

0.16285E+04 

0.0102 

N.O.F.:  Nodal  degree  of  freedom 
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It  may  be  noticed,  however,  that  boundary  nodes  behave  differenly  from  the 
interior  nodes  and  tend  to  have  greater  error  margins.  The  error  in  the  velocities 
of  the  boundary  nodes  is  due  to  the  complex  effects  of  friction,  temperature 
gradient,  and  changes  in  boundary  conditions.  In  spite  of  these  errors  we  still 
have  very  accurate  velocity  field  modeling  for  elements  within  the  boundary. 
Therefore,  in  the  present  design  scheme,  elements  at  the  boundary  are  generally 
not  selected  as  elements  of  interest. 

Based  on  the  above  example  and  other  similar  cases,  we  may  conclude  that 
the  state  space  model  is  capable  of  satisfying  the  design  requirements  and  can 
be  used  as  an  effective  numerical  model  for  design  of  the  process  parameters.  It 
was  also  observed  that,  overall,  the  state  space  model  was  found  to  give  more 
accurate  results  for  the  frictionless  case  as  compared  to  the  friction  case.  Also, 
the  results  of  the  isothermal  model  were  generally  found  to  be  better  than  those 
of  the  nonisothermal  state  space  model. 

6.2  Temperature  gradient  reduction 

One  of  the  primary  design  objectives  in  this  work  is  to  reduce  the  tempera¬ 
ture  gradient  at  the  boundaries,  because  this  would  alleviate  the  thermal  stresses 
induced  at  the  boundaries  due  to  the  die  chilling  effect.  This  example  illustrates 
the  influence  of  temperature  gradient  reduction  on  the  boundary  nodal  tem¬ 
peratures.  The  design  requirements  are  chosen  arbitrarily  just  to  illustrate  the 
method.  It  is  assumed  that  the  temperature  gradient  at  node  18  needs  to  be 
reduced,  and  the  effective  strain-rate  of  element  23  (£23)  is  to  be  maintained  at 
2.0  1/sec.  The  only  design  parameter  is  the  die  velocity.  The  process  is  started 
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with  a  relatively  low  die  velocity  of  0.5  in  /sec,  which  results  in  £23  around  0.25 
1  /sec.  This  gives  enough  room  to  change  the  die  velocity  by  large  amounts  to 
illustrate  the  effectiveness  of  the  current  method.  For  comparison  purposes,  two 
sets  of  designs  are  carried  out,  and  the  results  are  presented  in  Figs.  10  through 
12. 

The  first  design,  called  ‘Design  I’,  does  not  include  temperature  gradient 
reduction,  and  merely  involves  controlling  the  strain  rate.  The  output  variable 
of  Design  I  is  £23 •  The  second  design,  ‘Design  IF,  is  similar  to  Design  I,  but  in 
addition  involves  a  temperature  gradient  reduction  of  3%  at  every  time  step,  and 
has  node  18’s  temperature  added  into  the  output  vector.  Eqs.  (5.1)  and  (5.2)  are 
used  to  obtain  the  desired  £23  and  T\ s,  respectively.  The  weighting  matrices  are 
the  same  for  the  two  designs,  except  that  Design  II  has  an  additional  weighting 
factor  corresponding  to  the  temperature  output.  During  Design  I  die  velocity 
increases  at  a  steady  rate  till  £23  reaches  2.0  1/sec  at  a  stroke  of  0.23  inch.  In 
Design  II,  the  strain-rate  increases  faster  than  in  Design  I  because  of  the  effect 
of  temperature  gradient  reduction.  To  reduce  the  temperature  gradient  the  heat 
generation  rate  needs  to  be  faster  than  normal.  This  can  be  achieved  only  by 
increasing  the  die  velocity.  As  a  result,  the  strain  rate  also  increases  at  a  faster 
rate  and  reaches  the  required  value  quicker  than  in  the  previous  case.  Wffen 
the  stroke  approaches  0.13  inch,  a  sudden  jump  in  the  strain-rate  occurs  because 
the  workpiece  comes  in  contact  with  the  far  end  of  the  die.  Consequently,  £23 
overshoots  the  desired  value  of  2.0  1/sec.  The  die  velocity  is  now  reduced  to  bring 
l23  back  to  2.0  1/sec.  The  desired  value  is  finally  achieved  at  a  stroke  of  0.17 
inch.  Examining  the  temperature  profiles  in  Fig.  12,  it  is  observed  that  before 
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Effective  strain-rate  (1/sec 


Die  stroke  (*0.01in) 


Figure  10:  Temperature  gradient  reduction  -  Strain  rates 
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Die  velocity  (in/sec 


Die  stroke  (*0.01in) 


Figure  11:  Temperature  gradient  reduction  -  Die  velocities 
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a  stroke  of  0.20  inch,  the  temperature  gradient  shows  a  definite  reduction  in  the 
second  case.  But  after  that,  the  gradients  of  both  curves  are  nearly  the  same 
because  the  die  velocities  in  both  cases  are  almost  the  same.  At  the  end  of  the 
process,  the  temperature  of  node  18  is  1609° F  in  Design  I  and  1635 0 F  in  Design 
II,  a  26° F  increase  in  the  latter  case.  It  may  be  concluded  that  temperature 
gradient  reduction  results  in  an  increase  in  temperature  of  the  die-contacting 
boundary  nodes.  It  should  be  noted,  though,  that  this  increase  is  largely  a  result 
of  the  reduced  die-workpiece  contact  time  because  the  die  velocity  is  higher  in 
Design  II  than  in  Design  I.  This  example  shows  that  using  the  optimal  design 
scheme,  a  suitable  die  velocity  profile  may  be  generated  that  maintains  the  strain 
rate  of  critical  elements  at  the  desired  value  and  at  the  same  time  reduces  the 
temperature  gradient  at  the  die-workpiece  interface. 

6.3  Isothermal  design 


This  example  uses  the  same  engine  disk  forging  model  described  earlier.  An 
isothermal  design  is  presented  which  requires  the  effective  strain-rate  of  element 
23  to  be  maintained  at  0.1  1/sec.  The  starting  die  velocity  is  about  0.25  in/sec, 
corresponding  to  an  £23  value  of  0.125  1/sec.  Fig.  13  shows  that  this  starting 
£23  value  is  quickly  reduced  to  0.1  1/sec  because  the  design  scheme  calls  for  a 
reduction  in  the  die  velocity.  During  the  initial  stages  of  forging  the  change  in 
£23  is  relative  small,  but  around  a  stroke  of  0.13  inch  a  large  jump  (nearly  a  50% 
increase)  in  £23  may  be  noticed.  This  is  because  of  the  sudden  change  in  the 
deformation  mode  of  the  workpiece  as  it  comes  in  contact  with  the  far  end  of  the 
die.  However,  as  shown  in  Fig.  14,  the  design  scheme  is  robust  and  effective,  and 
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immediately  reduces  the  die  velocity  and  brings  the  effective  strain-rate  back  to 
0.1  1/sec  in  a  short  time  without  much  oscillation.  Near  the  end  of  the  process, 

,  though,  some  instability  is  observed,  causing  a  considerable  amount  of  oscillation 

in  both  the  die  velocity  and  strain-rate.  This  could  be  due  to  the  complex  nature 
of  forces  that  come  into  play  during  the  die  filling  process. 

* 

6.4  Temperature  range  reduction 

In  most  hot  die  forging  processes,  one  of  the  design  objectives  is  to  reduce 
the  temperature  range  in  the  billet  as  much  as  possible.  This  improves  the 
temperature  distribution  and  results  in  a  more  uniform  distribution  of  properties 
in  the  final  workpiece.  The  following  numerical  example  is  presented  to  illustrate 
how  the  above  objective  may  be  met  using  the  proposed  design  methodology. 
The  optimization  problem  is  stated  as: 

Minimize  :  | Tmaz  —  7min|  (®*^) 

Subject  to  :  Imax  <  2.5  1/ sec  (6*2) 

Design  variable :  Vj, 

where  Eq.  (6.2)  applies  a  constraint  that  requires  the  maximum  effective  strain- 
rate  in  the  billet  to  be  no  more  than  2.5  1/sec.  The  present  design  scheme  is 
used  to  generate  an  optimum  ram  velocity  profile  to  solve  the  problem  stated 
above.  The  optimal  results  are  then  compared  with  the  constant  velocity  results 
4  to  demonstrate  the  effectiveness  of  the  proposed  design  methodology. 

Through  finite  element  simulations  it  was  found  that  node  18  possesses  the 
lowest  temperature  in  the  billet  throughout  the  process.  This  is  because  node  18 
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remains  in  contact  with  the  die  from  beginning  to  end,  and  the  strain  rate  in  the 
vicinity  of  node  18  is  not  very  high.  The  temperature  range  in  Eq.  (6.1)  may  now 
be  approximated  to  ( Tmax  -  T18),  where  T18  is  the  temperature  of  node  18.  The 
design  scheme  has  the  ability  to  automatically  select  the  element  possessing  the 
largest  strain-rate  as  the  element  of  interest  and  design  an  optimal  die  velocity 
schedule  based  on  the  strain  rate  constraint  imposed.  Fig.  15  gives  the  strain- 
rate  schedule  for  the  element(s)  having  the  largest  strain  rate  in  the  workpiece 
at  any  given  time  step,  and  Fig.  16  gives  the  corresponding  optimal  die  velocity 
profile.  With  a  starting  value  of  2.1  in/sec,  the  designed  ram  velocity  gradually 
reduces  keeping  imax  around  2.5  1/sec  for  the  entire  process.  The  designed  ram 
velocity  profile  gives  the  largest  value  of  die  velocity  that  can  be  used  without 
violating  the  constraint  of  Eq.  (6.2).  Hence,  the  profile  also  gives  the  minimum 
contact  time  between  the  die  and  billet.  By  minimizing  the  contact  time,  the 
total  heat  loss  and  the  temperature  range  in  the  workpiece  are  also  reduced.  At 
the  end  of  the  process,  node  18’s  temperature  (the  lowest  temperature  m  the 
billet)  is  1574°  F  and  the  maximum  temperature  in  the  billet  is  1777°  F,  giving  a 
temperature  range  of  203 °  F. 

For  comparison,  the  above  example  was  also  simulated  using  a  constant  die 
velocity.  Because  tmax  gradually  increases  under  a  constant  die  velocity,  the 
largest  constant  die  velocity  which  can  be  used  without  exceeding  the  constraint 
is  1.25  in/sec.  By  examining  Figs.  15  and  16  it  may  be  observed  that  even 
though  for  most  of  the  process  lmar  is  far  below  2.5  1/sec  limit,  we  can  use  a 
constant  die  velocity  no  larger  than  1.25  in/sec;  otherwise  imax  would  violate  the 
strain  rate  constraint  towards  the  end  of  the  stroke. 


Figure  15:  Temp,  variance  reduction  -  Strain  rates 
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/sec) 


w 


* 


Figure  16:  Temp,  variance  reduction  -  Die  velocities 


94 


In  this  example  the  constant  velocity  is  lower  than  the  optimal  ram  velocity 
for  most  of  the  forging  process  (Fig.  16).  This  results  in  a  longer  die-workpiece 
contact  time  in  the  former  case.  At  the  end  of  the  simulation  it  may  be  observed 
that  node  18  has  a  temperature  of  1533 °F,  and  the  maximum  temperature  in 
the  billet  is  1775°  F.  The  temperature  range  thus  turns  out  to  be  242°  F.  Fig. 
17  clearly  illustrates  the  difference  in  temperature  profiles  (of  node  18)  for  the 
two  cases,  and  Fig.  18  gives  the  corresponding  load-stroke  curves.  There  is  only 
a  slight  increase  in  load  while  using  the  optimal  die  velocity  because  the  design 
scheme  is  based  on  the  minimum  energy  concept.  This  example  shows  that 
by  using  this  design  scheme,  an  optimal  die  velocity  profile  can  be  generated 
to  minimize  the  die-workpiece  contact  time  and  the  temperature  range  in  the 
deforming  workpiece.  Table  2  lists  some  of  the  performance  benefits  derived  by 
using  this  methodology. 


Table  2.  Performance  of  the  optimal  control  design  scheme 


Optimum  Velocity 

Constant  Velocity 

Max.  Temp.  (°F) 

1777 

1775 

Min.  Temp.  (°F) 

1574 

1533 

Temp.  Range  ( °F ) 

203 

242 

Range  Reduced  (%) 

16 

— 

Process  Time  (sec) 

0.300 

0.402 

Time  Reduced  (%) 

25.4 

— 
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6.5  Initial  die  temperature  design 


In  forging  processes,  selection  of  initial  die  temperature  is  an  important  but 
difficult  decision  to  make.  This  numerical  example  describes  how  the  proposed 
design  methodology  could  be  used  to  obtain  an  optimal  initial  die  temperature 
value.  The  design  problem  for  this  example  may  stated  as  follows: 


Objective  : 

Tig  >  1560° F 

(6.3) 

Subject  to  : 

I23  =  0.7  1/sec 

(6.4) 

Design  variables  :  Td  and,  Vd 


where  I23  is  the  effective  strain-rate  of  element  23  (critical  element  as  explained 
earlier),  and  Td  is  the  initial  die  temperature.  Because  node  18  (critical  node) 
possesses  the  lowest  temperature  in  the  billet,  Eq.  (6.3)  is  equivalent  to  raising 
the  temperature  of  the  whole  billet  above  1560°  F.  The  design  is  carried  out  in 
three  steps,  each  necessitating  an  independent  finite  element  simulation. 

Step  1.  Generating  the  desired  nodal  temperature  profile. 

In  this  step,  the  design  is  conducted  using  only  die  velocity  as  the  design  vari¬ 
able.  The  initial  die  temperature  used  is  600°  F.  The  purpose  of  this  step  is  to: 
(1)  check  whether  Eq.  (6.3)  can  be  satisfied  without  changing  the  initial  die  tem¬ 
perature,  and  (2)  generate  a  suitable  polynomial  representing  the  temperature 
path  of  node  18  that  satisfies  the  constraints  of  Eq.  (6.4).  The  design  is  started 
when  the  die  velocity,  Vd,  is  equal  to  1.4  in/sec,  corresponding  to  e23  equal  to 
0.698  1/sec.  Checking  the  results  at  the  end  of  the  simulation,  the  temperature 
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of  node  18  temperature  was  found  to  be  1519° F.  Figs.  20  through  22  give  the 
profiles  of  I23,  Vd  and  Tig,  respectively  for  this  example.  It  may  be  observed 
from  Fig.  20  that  the  strain  rate  constraint  of  Eq.  (6.4)  is  satisfied  for  the  entire 
process  under  the  optimal  die  velocity  (Vd)  in  Fig.  21.  At  the  same  time,  Fig. 
22  shows  that  the  temperature  constraint  of  Eq.  (6.3)  is  not  satisfied.  This  calls 
for  an  increase  in  the  initial  die  temperature  so  as  to  bring  the  temperature  of 
node  18  to  the  required  value.  To  represent  the  desired  temperature  profile,  a 
3rd  order  polynomial  is  used  to  fit  the  path  of  Tis  as, 

Tig  =  1735  -  641.374  +  355.5t2  -  25.lt3  (6.5) 

In  Eq.  (6.5)  the  first  order  term  is  the  dominant  term  influencing  the  temperature 
gradient.  To  raise  the  nodal  temperature,  therefore,  the  coefficient  of  the  first  or¬ 
der  term  has  to  be  changed  to  an  appropriate  value.  By  setting  Tig  =  1560°  F  at 
end  of  the  process,  the  new  first  order  coefficient  (from  Eq.  (6.5))  was  calculated 
as  -550.  Substituting  this  back  into  Eq.  (6.5)  we  have, 

Tig  =  1735  -  5504  +  355.5t2  -  25.lt3  (6.6) 

Eq.  (6.6)  gives  the  desired  temperature  profile  for  node  18  during  the  initial  die 
temperature  adjustment  design. 

Step  2.  Adjusting  the  initial  die  temperature 

With  Eq.  (6.6)  as  the  desired  temperature  profile,  the  design  in  step  1  is 
repeated  using  an  additional  design  variable  A Td,  which  is  the  initial  die  temper¬ 
ature  adjustment  parameter.  The  other  conditions  remain  the  same  as  in  Step 
1.  The  value  of  A  Td  at  every  time  step  is  calculated  using  the  proposed  design 
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scheme  and  presented  in  Fig.  19.  By  using  Eq.  (5.8)  to  average  the  A Td  values 
at  each  simulation  time  step,  a  mean  initial  die  temperature  adjustment  value  of 
172° F  is  obtained  (represented  by  the  dashed  line  in  Fig.  19).  The  actual  initial 
die  temperature  which  satisfies  all  the  design  constraints  is  thus  112°  F. 

Step  3.  Redesign  using  the  new  initial  die  temperature 

This  step  is  merely  a  repeat  of  step  1  with  the  new  initial  die  temperature 
value  of  112°F  (instead  of  600°.F).  From  Figs.  20  and  21,  the  profiles  of  the 
die  velocity  and  £23  are  similar  before  and  after  initial  die  temperature  design. 
This  shows  that  changing  the  initial  die  temperature  does  not  affect  the  strain- 
rate  of  the  control  element  significantly.  The  temperature  of  node  18,  though, 
is  very  sensitive  to  the  change  in  initial  die  temperature  as  observed  in  Fig.  22. 
With  a  112°  F  higher  initial  die  temperature,  the  final  temperature  of  node  18 
is  1557° F,  very  close  to  the  design  requirement  of  1560°F  (a  difference  of  about 
7.5%).  The  temperature  gradient  of  the  boundary  nodes  (node  18  in  particular) 
is  thus  reduced.  At  the  same  time,  it  was  also  observed  that  the  temperature 
range  was  reduced  from  256°  F  to  218°F,  a  drop  of  about  14%. 

This  example  shows  that  the  proposed  method  for  initial  die  temperature 
(adjustment)  design  works  well,  and  can  be  effectively  used  to  reduce  the  tem¬ 
perature  range/temperature  gradient  in  the  deforming  workpiece.  Noting  that 
the  increase  in  boundary  nodal  temperatures  is  not  linearly  proportional  to  the 
increase  in  initial  die  temperature,  solving  a  design  problem  such  as  the  one 
posed  in  this  example  would  be  very  difficult  and  cumbersome  without  a  sys¬ 
tematic  methodology  and/or  strategy  like  the  one  described  in  this  section.  One 
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of  the  major  advantages  of  the  method  described  here  is  that  the  design  can  be 
performed  in  three  simulation  steps  instead  of  going  through  a  number  of  trial 
and  error  steps  as  is  done  normally.  This  could  save  a  considerable  amount  of 
time,  effort,  and  money  during  the  entire  design  process. 


4 
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Figure  19:  Initial  die  temperature  adjustment  * 
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A  Figure  22:  Initial  die  temp,  design  -  Node  18  temp,  profiles 
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CHAPTER  7 


Reduced  Order  State  Space  Models 


* 


7.1  Introduction  to  model  reduction 

In  general,  large  scale  finite  element  models  are  used  for  simulating  and 
analyzing  metal  forming  processes  to  capture  in  detail  information  regarding 
the  thermo-mechanical  behavior  of  the  deforming  material.  The  condensation 
scheme  explained  earlier  reduces  and  transforms  the  finite  element  model  to  a 
state  space  model  by  retaining  the  tangential  nodal  velocities  and  temperatures 
of  the  die  contacting  nodes,  and  the  nodal  degrees  of  freedom  of  the  critical 
element  (in  the  deforming  material)  as  states.  Despite  the  condensation  process, 
a  large  model  results  when  analyzing  large  scale  practical  forming  problems. 

Further,  for  on-line  control  applications,  it  is  desirable  to  have  smaller  models 
capable  of  accurately  representing  the  thermo-mechanical  behavior  of  the  metal 
forming  process.  Since  the  state  space  model  is  built  using  the  converged  finite 
element  solution  at  each  time  step,  it  may  be  considered  a  reasonably  accurate 
representation  of  the  metal  forming  process  [66,67]. 

The  highly  nonlinear  finite  element  model  of  the  metal  forming  process  is  f 

linearized  and  modeled  in  the  state  space  form  as, 

x  =  Ax  +  Bn 
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y  —  Cx 


(7.1) 


where  A  is  the  plant  matrix,  B  is  the  input  matrix,  and  C  is  the  output  matrix  to 
^  convert  the  state  vector  x  into  the  desired  output  y .  u  is  the  input  vector,  which 

in  this  case  is  the  ram  velocity.  Eq.  (7.1)  is  derived  from  the  stiffness  matrix 
4  and  force  vectors  of  the  metal  forming  process.  The  state  space  models  are 

built  using  the  converged  solution  at  the  end  of  each  simulation  time  step.  The 
conversion  of  the  finite  element  equations  to  state  space  form  takes  into  account 
the  symmetric  boundary  conditions  and  the  other  relevant  boundary  conditions 
of  the  physical  system.  The  force  vector  used  includes  the  effects  due  to  friction 
and  reaction  forces.  The  reaction  forces  are  due  to  the  symmetric  boundary 
conditions  applied  during  the  finite  element  modeling.  The  nodal  velocities  of 
the  boundary  nodes  and  the  nodal  velocities  of  the  element  of  interest  (whose 
strain  rate  is  to  be  controlled)  are  considered  as  states.  The  elemental  strain  rate 
is  the  output,  and  it  is  maintained  within  a  predetermined  and  prespecifed  range. 
The  output  matrix,  C,  is  defined  as  the  matrix  to  convert  the  nodal  velocities 
into  the  strain  rate  of  the  element  of  interest.  The  control  problem  is  posed  as 
a  tracking  problem,  and  since  the  output  is  to  be  tracked,  the  objective  of  the 
control  law  is  to  minimize  the  performance  index  in  such  a  way  that  the  error  is 
minimized  and  the  end  conditions  are  met. 

The  aim  of  reducing  the  original  model  is  to  get  an  equivalent  linear  state  s- 
pace  system  of  a  dimension  considerably  smaller  than  that  of  the  original  system. 
*  Since  the  control  problem  is  a  tracking  problem,  the  reduced  model  should  nec¬ 

essarily  be  controllable.  In  addition,  the  reduced  system  should  be  observable 
4  in  order  to  implement  this  methodology  as  a  feedback  on-line  control  scheme. 
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The  methodology  for  obtaining  reduced  order  state  space  models  and  their  sub¬ 
sequent  use  in  designing  optimal  die  velocity  schedules  is  decribed  in  the  flow 
chart  of  Fig.  23. 

Reduced  order  models  can  generally  be  obtained  using  several  methods  avail¬ 
able  in  relevant  literature.  The  following  are  two  broad  classifications  of  currently 
available  order  reduction  techniques: 

i.  The  first  classification  uses  the  mode  concept  of  linear  models.  The  aim  here 
is  to  keep  the  dominant  eigenvalues  of  the  full  model,  i.e.,  the  eigenvalues  most 
influential  on  the  system  dynamic  behavior.  The  modal  analysis  method  and 
the  aggregation  method  belong  to  this  classification. 

ii.  The  second  classification  includes  the  approaches  based  on  preserving  the 
same  degree  of  controllability  and  observability  product  (oy)  as  in  the  full 
model.  The  Balanced  Model  Reduction  (BMR)  method  belongs  to  this  cate¬ 
gory,  where  the  reduction  is  applied  by  deleting  the  smallest  singular  values, 
o-{,  of  the  balanced  state  space  configuration. 


» 
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7.2  Aggregation  method 


Aggregation  techniques  replace  the  initial  system  model  with  a  reduced  order 
model  that  preserves  the  dominant  system  characteristics.  In  this  method,  the 
important  assumption  is  that  the  reduced  state  vector  (xT)  and  the  full  state 
vector  x  are  related  linearly  by  an  aggregation  matrix,  P,  of  dimension  (mxp), 
having  rank  tti.  Here,  ttl  is  the  number  of  states  in  the  reduced  order  model,  and 
p  represents  the  number  of  states  in  the  full-size  model.  The  aggregated  model 
is  not  a  part  of  the  full  state  model.  Instead,  it  is  a  linear  combination  of  the 
states  of  the  full  model  [55].  Some  salient  features  of  this  technique  are  described 
below.  In  this  method  it  is  assumed  there  exists  an  aggregation  matrix  (P)  such 
that, 

xr  =  Px  (7-2) 

where  xr  is  the  reduced  state  vector,  and  x  is  the  full  state  vector.  Also,  xe 
xTe  Rm  and  Pe  R(mxp).  Now,  the  aggregated  state  space  model  can  be  con¬ 
structed  as, 

xr  =  ATxr  4-  Bru 

yr  =  Crx7  (7.3) 

where  Are  R(mxm\  BTeR<mxl)  and  Cre  R(lxm)  are  the  aggregated  plant,  input 
and  output  matrices,  respectively.  Since  the  input  and  output  are  scalars  (for  the 
isothermal  case),  Br  and  Cr  take  the  sizes  described  above.  In  order  to  derive 
the  aggregated  model  from  the  full  state  model,  the  relationship  between  the  full 
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model,  aggregated  model,  and  the  transformation  matrix  must  be  determined. 
By  substituting  Eq.  (7.2)  in  Eq.  (7.3),  we  have, 

Px  =  ArPx  +  Bru  (7.4) 

By  premultiplying  Equation  (7.1)  with  P  and  equating  the  coefficient  matrix  of 
x  and  the  coefficient  matrix  of  u,  with  those  of  Equation  (7.4),  we  get, 


* 


4- 


ArP  =  PA 

Br  =  PB  (7.5) 

The  selection  of  P  is  important  to  establish  a  proper  relationship  between  the 
full  state  model  and  aggregated  model.  Generally,  this  matrix  is  computed  based 
on  the  contribution  of  the  eigenvalues  to  the  output  vector.  If  A  represents  the 
eigenvalues  of  A,  and  v  represents  their  associated  eigenvectors,  then, 

Av  =  Av  (7.6) 

Premultiplying  both  sides  of  Equation  (7.6)  by  the  aggregation  matrix,  P, 

P  ( Av)  =  APv  (7.7) 

Substituting  for  PA  from  Equation  (7.5)  into  Equation  (7.7), 

Ar(Pv)  =  A(Pv)  (7.8) 

This  equation  shows  that  the  eigenvalues  of  the  matrix  Ar  are  m  aggregated 
eigenvalues  of  the  matrix  A.  The  same  is  true  in  case  of  the  eigenvectors  also. 
Hence,  it  may  be  seen  that  the  selection  of  the  aggregation  matrix  significantly 
influences  the  performance  of  the  reduced  model  with  respect  to  the  full  model. 
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Since  the  aggregation  matrix  is  a  transformation  matrix  which  is  selected  based 
on  eigenvalues,  the  selection  is  carried  out  using  the  modal  matrix  of  the  full 
system.  The  canonical  form  of  the  states  is  expressed  as, 

x  =  Mw  (7.9) 


where  w  is  the  transformed  state  variable  in  eigenspace  and  M  is  the  modal 
matrix  containing  the  eigenvectors  of  the  matrix  A.  By  substituting  Equation 
(7.9)  into  Equation  (7.1),  the  full-size  state  space  system  in  the  transformed  form 
is  written  as, 


where, 


w  =  Aw  +  Tu 

(7.10) 

=  M"1AM=(Ao1  a°2) 

(7.11) 

r  =  m_1b  =  (p’) 

(7.12) 

Here,  A  is  a  diagonal  matrix  of  the  eigenvalues  of  A,  Ai  contains  the  dominant 
eigenvalues  and  A2  is  made  up  of  the  non-dominant  eigenvalues  of  the  system. 
Because  of  the  properties  of  the  stiffness  matrix,  the  plant  always  has  distinct 
eigenvalues.  TieR(mxm)  and  T2e  R(P-m)x(P”m)  are  the  transformed  forms  of  the 
input  matrix  of  the  full  state  model.  In  the  partitioned  form  the  full  state  space 
system  is  written  as, 

'  Ai  0  WwE 


(£M 


0  A2J  l.w2)  +  (r2)“ 


(7.13) 


From  this,  the  decomposition  of  the  transformed  states  into  wi  reduced  states  is 
performed  as  described  below: 


wi  =  A]  wi  -f  Tiu 


(7.14) 


where  Wi  represents  states  that  are  to  be  retained  for  the  reduced  model  of  size 
m,  and  is  defined  as  wj  =  (Im  0  )w.  Also,  substituting  for  w  from  Equation 
(7.9),  we  get  Px  =  (Im  0  )M-1x.  Now,  the  transformation  matrix  P  can  be 
defined  as, 

P  =  (Im  0  )M-1  (7.15) 

Because  the  goal  of  model  reduction  is  to  reduce  the  size  of  the  original  model, 
but  still  obtain  results  comparable  to  that  of  the  full  model,  the  output  matrix 
Cr  also  must  be  defined  in  terms  of  the  aggregation  matrix.  Hence,  the  reduced 
output  matrix  is  given  by, 

Cr  =  CP+  (7.16) 

where  P+  is  the  pseudo  inverse  of  the  aggregation  matrix  P. 


7.3  Modal  analysis  method 

In  the  aggregation  method  the  states  (reduced)  are  a  linear  combination  of  the 
full  state  model.  The  modal  analysis  method,  unlike  the  aggregation  method,  is 
based  on  transforming  the  system  matrices  in  such  a  way  that  the  characteristics 
of  the  entire  system  (full-size)  are  included  in  the  reduced  model.  The  full  state 
model  is  represented  in  eigenspace  using  the  relation  in  Equation  (7.9),  where 
the  transformation  is  done  using  the  modal  matrix,  M.  Equation  (7.1)  is  valid 
for  one  simulation  time  step  of  the  metal  forming  process.  Because  the  system 
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(Eq.  (7.1))  is  now  linear  and  time  invariant  for  that  time  step  (interval),  it  can 
be  partitioned  in  the  following  form, 

(*;)=(£  (7-i7) 
where  xieRm  is  the  state  vector  whose  modes  are  to  be  retained,  x2eR^_m^  is 
the  state  vector  whose  modes  are  to  be  neglected,  and  m  is  known  a  •priori.  The 
matrices  A1}  A2,  A3,  A4,  Bx  and  B2  are  constant.  Accordingly,  the  modal 
matrix,  M,  can  also  be  represented  in  the  partitioned  form  as, 


M  = 


/Mi 

vm3 


m2\ 

M  J 


(7.18) 


From  Eq.  (7.9),  x  =  [xTi  xT2]T  is  the  original  form  and  w  =  [wTi  wT2]T 
is  the  transformed  form.  In  the  partitioned  form,  wa  includes  the  states  of 
dominant  eigenvalues  and  w2  includes  states  of  non-dominant  eigenvalues  of  the 
system.  Now,  the  original  system  is  represented  in  the  transformed  coordinates 
as  explained  in  equations  (7.10)  -  (7.12).  The  states  in  eigenspace  are  further 
expanded  using  the  matrix  properties  and  the  relationship  between  the  matrices 
of  M  and  their  inverses.  The  dominant  and  non-dominant  states  are  expressed 
as, 


wi  =  Mi  1  [AjMj  -+-  A2M3]  wi 

+  (Mi  -  M2M4-1M3)  Bi  -  Mf'Mj  (M4  -  M3Mi"1M2)  '  B2  u  (7.19) 

w2  =  M4-1  [A3M2  +  A4NI4]  w2 

+  f-M4-1M3(M4  -M2M4-1M3)_1Bi  +  (M4-M3Mf1M2)  1  B2  u  (7.20) 


From  the  above  two  expressions  the  reduced  order  models  may  be  derived.  Based 
on  the  modal  analysis  technique,  there  are  three  methods  by  which  model  reduc¬ 
tion  is  done: 


JT 


* 


* 
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1.  Davison’s  method 


■* 


A 


2.  Marshall’s  method 

3.  Nicholson’s  method 

All  the  above  methods  retain  the  dominant  eigenvalues  of  the  original  system. 
Further,  all  three  methods  assume  that  all  the  non-dominant  eigenvalues  are  in 
A2.  Because  the  design  is  a  tracking  problem,  the  output  matrix  is  also  modified 
accordingly,  to  get  reasonable  results.  Some  distinct  features  of  these  methods 
are  discussed  below. 


7.3.1  Davison’s  method 

This  method  makes  use  of  the  fact  that  A2  contains  only  non-dominant 
eigenvalues.  w2  is  derived  from  the  approximation  w2  =  0  and  is  replaced  in 
terms  of  Xi  and  u.  By  substituting  this  approximation  back  into  Equation  (7.10) 
and  rearranging  the  terms,  the  reduced  system  is  written  as, 

At  =  Ai  +  A2M3Mi_1 

Br  =  Mi  (Mi  -  M2M4_1M3)  1  Bi  -  M2  (M4  -  M3Mi_1M2)  B2  (7.21) 

f  Because  the  output  of  the  reduced  model  depends  on  the  retained  states  xi,  the 

output  matrix  is  also  partitioned  as, 

*  y  =  [Ci  c2,  (*‘)  =  [Ci  c2]  (mJ  EXS) 
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where  Ci  and  C2  are  the  partitioned  output  matrices  spanning  the  dominant 
and  non-dominant  states,  respectively.  By  rearranging  and  substituting  for  the 
non-dominant  states,  X2,  the  output  is  written  as, 


y  =  [Ci  +  C2M3M11]  xi  +  {C2  -  M4]  A2_1r2}u  =  Crxi  +  Dru 

(7.23) 


7.3.2.  Marshall’s  method 

This  method  is  derived  from  Equation  (7.17)  and  equations  (7.10)  -  (7.12). 
From  Equation  (7.17),  Xi  can  be  defined  as, 

xi  =  Ajxi  +  A2X2  -I-  Biu  (7-24) 

By  partitioning  Equation  (7.9)  and  extracting  the  terms  related  to  X2,  we  get, 

X2  =  M3W1  +  M4W2  (7.25) 

Now,  rearranging  Equation  (7.25),  the  dominant  states  in  the  modal  coordinates 
are  obtained  as, 

wi  =  Mi"1xi-MrIM2W2  (7.26) 

Using  the  condition  W2  =  0,  W2  can  be  determined  in  the  modal  coordinates  as, 

W2  -  -A2-1r2U  (7.27) 

By  substituting  Equations  (7.26)  and  (7.27)  in  Equation  (7.25),  and  substituting 
the  result  in  Equation  (7.24),  the  reduced  system  is  derived  as, 


A  r  =  Ai  +  A2M3M1  1 


> 


M. 


Br  =  Bi  -  A2  [M4  -  M3Mf1M2]  A2_1r2  (7.28) 

In  modal  coordinates  the  output  is  written  as  in  Equation  (7.22).  Now,  substi¬ 
tuting  for  w2  from  Equation  (7.27),  the  output  (in  terms  of  reduced  states)  may 
be  written  as, 

y  =  [CiMi  +  C2Ms]  wi  -  [CiM2  +  C2M4]  A2_1r2u  =  Crwx  +  D Tu  (7.29) 


7.3.3.  Nicholson’s  method 

In  this  method  the  reduced  system  is  derived  by  completely  neglecting  the 
effect  of  w2  on  the  system.  By  assuming  that  the  non-dominant  states  do  not 
have  any  influence  on  the  dominant  states,  and  substituting  for  x2  in  terms  of 
Xl,  the  reduced  model  is  constructed.  Based  on  this  assumption,  Equation  (7.18) 
is  modified  as, 

(SHE  °.)GD  (7-30) 

By  substituting  Equation  (7.30)  in  Equation  (7.17)  and  rearranging,  we  get, 

xi  =  MiAiM^1xi  +  Mil'll*  (7-31) 

The  reduced  model  may  now  be  described  as, 

A  7*  =  Ai  +  A2MsMi  1 

Br  =  Mi  ( Im  0mx(p_m) )  M_1B  (7.32) 
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The  neglected  (non-dominant)states  are  given  by, 

x2  =  MjM^xi 


(7.33) 


The  output  matrix  for  the  reduced  system  may  be  derived  as, 


y  =  [Ci  -f  Xi  —  Crxi 


(7.34) 
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7.4  Balanced  model  reduction  (BMR)  method 

7.4.1  Construction  of  the  balanced  model 

The  BMR  method  belongs  to  the  second  classification  of  order  reduction 
methods  described  in  section  7.1.  Generally,  it  is  desirable  to  have  a  system 
that  is  both  controllable  and  observable,  called  the  balanced  system.  The  BM- 
R  method  is  a  state  selection  procedure  in  the  grammian-balanced  coordinate 
system  of  the  full  state  space  model.  The  grammian-balanced  system  is  called 
the  internally  balanced  system  of  the  full  system.  The  reduced  model  is  derived 
by  retaining  only  the  most  controllable  and  observable  states  of  the  internally 
balanced  system.  The  internally  balanced  system  is  constructed  using  a  trans¬ 
formation  matrix,  [Tb],  which  exists  only  if  the  plant  matrix  of  the  full  state 
space  model  is  asymptotically  stable.  The  state  space  model  for  metal  forming 
processes  is  always  stable  because  the  plant  matrix  is  negative  definite,  as  may 
be  observed  from  chapter  2.  If  the  transformation  matrix  exists,  then,  the  states 
in  the  balanced  coordinate  system  are  given  by, 

xb  =  Tb-1x  (7.35) 

Now,  the  balanced  state  space  form  of  Equation  (7.1)  is  expressed  as, 

xb  =  Abxb  +  Bbu  +  Wb 


y  =  Cbxb 


(7.36) 


where  Ab  =  Tb1ATb,  Bb  =  Tb  X,  Cb  —  CTb  and  Wb  —  Tb  W.  A 
system  is  called  internally  balanced  or  grammian  balanced  if, 

Pc  =  Po  =  s  =  diag  [<ri2,  £r22,  •  •  •  ,o-n2]  (7-37) 

where  the  cq’s  are  Hankel  singular  values  of  the  balanced  system  and  E  is  a  diag¬ 
onal  matrix  containing  <tj2.  Pc  and  P0  are  the  controllability  and  observability 
grammians  of  the  system,  which  are  the  solutions  of  the  following  Lyapunov 
equations: 

AbPc  +  PcATb  +  BbBTb  =  0 

and 

ATbP0  +  PoAb  +  CTbCb  =  0  (7.38) 

The  balanced  state  space  model  (Equation  7.36)  is  actually  the  controllable  sub¬ 
space  of  the  condensed  (full-size)  state  space  model  of  the  metal  forming  process 
(Equation  7.1). 

7.4.2  Construction  of  the  transformation  matrix  [57] 

The  controllability  and  observability  grammians  (Pc  and  PQ)  of  the  full  state 
space  model  may  be  obtained  as  solutions  of  the  following  Lyapunov  equations: 

APC  +  PcAt  +  BBt  =  0 

and 

AtP0  +  PoA  +  CTC  =  0.  (7.39) 
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By  factorizing  the  controllability  grammian  Pc  using  singular  value  decomposi¬ 
tion,  we  have, 

Pc  =  UlSiUi1  (7.40) 

From  the  properties  of  singular  value  decomposition,  the  product  U1U1  results 
in  an  identity  matrix,  and  Si  is  a  diagonal  matrix  having  the  singular  values  of 
the  matrix  Pc  as  its  diagonal  terms.  Si  will  normally  have  nc  non-zero  values. 
From  the  non-zero  values  of  Si,  the  controllable  subspace  of  the  system  can 
be  extracted  by  partitioning  the  matrix  U  as  [Uu  U12],  where  UneRn  x  nc 
spans  the  controllable  subspace  and  Ui2eRn  x  nc^  spans  the  uncontrollable 
subspace  of  the  system.  If  all  the  states  are  controllable,  then  nc  =  n.  Using  the 
information  specified  above,  the  controllable  subspace  may  be  retained  as, 

Ti  =  UiiSu*  (7-41) 

By  factorizing  the  product  TiTP0Ti  using  singular  value  decomposition 
[56,57,68], 

TiTP0Ti  =  U2E2U2T  (7.42) 

where  the  non-zero  portion  of  the  singular  value  matrix,  S2 ,  spans  the  observable 
part  of  the  controllable  subspace  of  the  system,  represented  as  X22.  By  partition¬ 
ing  U2  similar  to  the  operation  in  Equation  (7.40),  the  controllable-observable 
subspace  is  constructed  as, 

T2  =  U2iS22-5  (7.43) 

From  Ti  and  T2,  the  transformation  matrix  for  balanced  model  reduction  (Tb) 
is  derived  as, 


Tb  =  TiT2 


(7.44) 


7.4.3  Construction  of  the  reduced  model 


As  explained  by  Moore  [68],  in  cases  where  a2m  is  considerably  larger  than 
°r2{m+l))  the  subspace  xT  behaves  like  it  is  both  controllable  and  observable.  If 
the  mechanics  of  Kalman’s  minimal  realization  theory  are  applied  to  the  inter¬ 
nally  balanced  model,  with  xr  used  as  a  working  approximation  of  x then,  the 
resulting  lower  order  model  is  generically  stable  and  internally  balanced.  Now,  by 
truncating  the  smallest  singular  values  (<7;)  of  the  balanced  system,  the  reduced 
model  may  be  defined  by, 

Ar  =  (  Im 

Br  =  (  Im  0  )  Bfc 

c,  =  Ct  (I0m)  (7.45) 

7.5  Control  law  for  reduced  order  models 

The  performance  of  reduced  order  models  can  be  measured  using  many  cri¬ 
teria.  The  objective  of  model  reduction  in  this  case  is  to  design  optimal  input 
velocities  while  using  a  smaller  state  space  system.  Hence,  comparing  the  die 
velocities  designed  using  the  full  and  reduced  models  is  a  good  measure  to  e- 
valuate  the  performance  of  the  reduced  order  models.  The  control  law  used  for 
the  full-size  state  space  model  [66,67]  cannot  be  directly  applied  to  the  reduced 
models  (Equations  (3.23)  and  (3.29)),  because  the  reduced  models  have  an  extra  > 
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term  (D)  in  the  output  equation.  In  this  section,  a  control  law  suitable  for  use 
with  reduced  order  models  is  developed.  The  input  velocities  are  designed  using 
a  finite  time  control  law  with  output  tracking. 

For  a  time-invariant  linear  system,  the  state  space  representation  may  be 
given  by, 

x  =  f  [x(t),u(t),t]  =  Ax  +  Bu  +  W  and  y  =  Cx  +  Du  (7.46) 

where  A,B,'W,C,D,x  describe  the  full  system  and  Ar,Br,1W^r5Dr,Dr,xr  de¬ 
scribe  the  reduced  system.  The  design  is  done  off-line  as  explained  earlier.  The 
problem  is  posed  as  an  output  tracking  problem  and  the  output  tracked  is  the 
strain  rate  of  the  element  of  interest,  td  is  the  desired  strain  rate  trajectory  to 
be  tracked,  and  y  is  the  output  of  the  system.  The  error  term  (e)  is  given  by 
e  =  id  —  Cx  —  Du.  This  error  signal  is  used  as  the  state  feedback  while  designing 
the  velocity  schedule.  For  control,  it  is  desirable  to  bring  the  system  from  a  given 
initial  state  xQ  to  an  acceptable  terminal  state  by  designing  a  suitable  control 
input  u.  The  control  is  achieved  by  minimizing  a  performance  index  consisting  of 
positive  definite  quadratic  terms  for  the  terminal  condition,  tracking  error,  and 
control  input,  respectively.  This  is  represented  mathematically  as, 

tf 

min  J  =  4>  [x  (tf )  ,  tf]  +  J  L  [x(t)  ,u  (t)  ,t]  dt  (7-47) 

to 

where  ^>[x(tf),tf]  =  |eT(tf)Se(tf)  and  L[x(t),u(t),t]  =  |[eT(t)Qe(t)  +  u  Ru]. 
S  and  R  are  positive  definite  weighting  matrices  which  weight  the  terminal  condi¬ 
tion  and  the  control,  respectively,  and  Q  is  a  positive  semi-definite  matrix  which 
weights  the  states.  Here,  the  problem  is  to  design  the  control  input,  u(t),  such 


that  the  Equation  (7.47)  is  minimized.  To  find  an  optimal  u(t),  the  system  dif¬ 
ferential  equation  (Eq.  (7.46))  is  combined  with  Equation  (7.47)  by  means  of  a 
Lagrange  multiplier  function,  A(t),  as  shown  below: 
tf 

J  =  <£[x(tf),tf]  +  |  [L [x(t),u(t),t]  +  AT(t){f  [x(t),u(t),t] -x}]  dt 
to 

(7.48) 

A(t)  is  a  function  of  the  states  x(t).  For  numerical  convenience  a  scalar  function 
H  called  the  Hamiltonian,  is  defined  as  [69]: 

H [x  (t) ,  u  (t) ,  t]  =  L  [x  (t) ,  u  (t) ,  t]  +  AT  (t)  f  [x  (t) ,  u  (t) ,  t]  (7.49) 

By  substituting  Equation  (7.49)  in  Equation  (7.48)  and  integrating  the  resulting 
equation,  the  function  to  be  minimized  is: 

tf 

J  =  0(x(tr),tf]  -  A‘  (tf)x(tf)  +  V (to)x(t„)  +  J {fl(x(t),u(t),t]  +  A*(t)x(t)}dt 

(7.50) 


k 


Now  consider  the  variation  of  J  due  to  variations  in  the  control  u(t),  for  the  fixed 
time  interval  tQ  to  tf.  Mathematically, 

tr 

—  /*nrl  4-  [AT^x1  4- 

l\0X 


S J  = 


+  [AT<X],„  +  /  t 


dH 


dx 


+  A‘ 


/3H 

Sx  H — jr— <5u}dt  (7.51) 

ax 


to 

By  setting  the  coefficients  of  Sx  to  zero  in  the  above  equation,  the  necessary  and 
sufficient  conditions  for  minimizing  the  performance  index  can  be  derived  as, 

dH 


AT(t)  = 


dx 


dH  n 
and  — —  =  0 

au 


(7.52) 


The  corresponding  boundary  condition  may  be  given  as, 

AT  ft  1  ^ 


(7.53) 
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The  optimal  u(t)  which  minimizes  the  performance  index  (Eq.  (7.47))  can  be 
obtained  by  solving  Equations  (7.46),  (7.52),  and  (7.53)  simultaneously.  The 
Hamiltonian  function,  H ,  is  now  derived  as, 

H  =  l  jeTQe  +  uTRu]  +  AT  [Ax  +  Bu  +  W] 

=  i  jZQZ  -  xtCtQZ  -  utDtQZ  -  ZQCx  +  xTCTQCx  +  uTDTQCx] 

+-  [-ZQDu  +  xtCtQDu  +  utDtQDu  +  uTRu]  +  AT  Ax  +  ATBu  +  ATW 

(7.54) 

where  Z  is  the  desired  strain  rate  td.  Taking  the  partial  derivative  of  Equation 
(7.54)  with  respect  to  x  and  setting  the  result  equal  to  -A,  one  of  the  Euler- 
Lagrange  equations  may  be  constructed  as, 

A  =  -AtA  -  CtQCx  +  CtQZ  -  CtQDu  (7.55) 

Now,  taking  the  partial  derivative  of  Equation  (7.54)  with  respect  to  u  and 
setting  the  result  equal  to  zero,  the  optimal  control  law  may  be  derived  as, 

u(t)  =  -  [DTQD  +  r]_1  [DTQZ-DTQCx(t)  +  BTA(t)]  (7.56) 

Substituting  for  u(t)  from  Equation  (7.56)  into  Equation  (7.46),  we  have, 

x(t)  =  Hix(t)  +  HjA(t)  +  Vx(t)  (7-57) 

where, 

Hi  =  A  +  B  [dtQD  +  r]_1DtQC 

H2  =  -B  [DTQD  +  R]"1  Bt 

Vi  (t)  =  W-B  [DTQD  +  R]-1DTQZ(t) 
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Substituting  for  u(t)  from  Equation  (7.56)  into  Equation  (7.55),  we  get, 


A  (t)  =  H3x(t)  +  H4A(t)  +  V2(t)  (7.58) 


where, 

h3  =  -ctqc  -ctqd  [dtqd  +  r]  dtqc 
h4  =  -at  +  ctqd  [dtqd  +  r]'1  bt 

V2  =  CTQZ  +  CTQD  [DTQD  +  R]_1DTQZ(t) 
Combining  Equations  (7.57)and  (7.58),  we  get, 


/x(t)\  /Hi 
\  A  (t) )  \H3 


h2 

h4 


+  (Vim 


(7.59) 


There  are  2n  number  of  equations  in  Equation  (7.59),  where  n  is  the  number  of 
states  in  Equation  (7.46).  Now,  the  solution  of  Equation  (7.59)  may  be  obtained 


as. 


(m$)  “  *(*'-*>(*$)  +  *(tr,t)jW.t)(£$)dr  (7.60) 

to 

where  is  the  state  transition  matrix  (explained  earlier).  For  convenience, 
Equation  (7.60)  can  be  rewritten  as, 


(7.61) 


The  two  terms  on  the  right  side  of  Equation  (7.61)  are  computed  using  the 


Pade  approximation  technique  [63].  In  Equation  (7.61)  there  are  two  unknowns 
(x(tf )  and  A(t))  out  of  which  only  A(t)  is  needed  for  the  design  of  the  system 
input  (Equation  7.56).  To  solve  for  these  unknowns  2n  boundary  conditions  are 

needed.  These  are  specified  as, 

at  t  =  t0,  x(tD)  is  known 


at  t  =  tf ,  A  (tf) 


d<f> 1 

3x 


t 


k 
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By  partitioning  *  in  Equation  (7.61)  into  four  parts,  Equation  (7.61)  can  be 
rewritten  as, 

/  x  ( tf )  \  ($11  *12  \ /x(t)\  (giW'l  (7.62) 

(am)  =  ($21  $22  )  (  A(t)  /  +  V  g2  (t )  / 

From  the  boundary  condition  (Equation  7.53)  and  the  definition  of  d,  A(tf)  may 

be  defined  as: 

d 

A(tf)= 


{1  [zT  -  xTCT  -  uTDt]  S  [Z  -  Cx  -  Du)} 


5x(tf) ’"2 
_CT  _ 

ax 


S  [Z  -  Cx  -  Du]  +  [ZT  -  xTCT  -  utDt]  S 


1 

du 

Is 

L  J 

=  — C*SZ  +  CTSCx  +  CtSDu  +  [DtSDu  +  DTSCx  DTSZ] 


(7.63) 


Using  the  definition  of  u  from  Equation  (7.56), 

=  CTQDk~T  =  a 


3ut 


dx 


(7.64) 


where  k  =  DtQD  +  R.  Substituting  Equation  (7.64)  in  Equation  (7.63)  and 
grouping  the  terms,  we  have, 

A  (tt)=  [— CTS  —  aDTs]  Z  +  [CTSC  +  aDTSC]  x  +  [cTSD  +  aDTSD]  u 

=  nZ  +  r2X  +  r3u  (7  65) 

Again  substituting  for  u  from  Equation  (7.56),  we  have, 

A  (tf)  =  nz  +  T2X  +  r,  [— k-1DtQCZ  +  x-1DTQCx  -  (t)]  (7.66) 

Rearranging  and  regrouping  the  above  equation,  we  get, 

[i  +  r.rT'B1]  A(tj)  =  [n  -  t,k-1DtQC]  Z  (tf)  +  [n  +  r3x-1DTQC]  x(tf) 

(7.67) 


127 


or, 


A(tf)=[l  +  T3/c  1BTj  {[ri-r3/c  1DtQ]  Z  (tf)  +  [r2  +  t3k  1DtQC  x(tf)} 
=  T7xZ  (tf)  +J?2x(tf) 


where 


771  =  [l  +  r3/c  1BT]  [ti  -  t3k  1DTQ 
772=  [l  +  r3/c_1BT]  [t2  +  t3k'1DTQC 


By  expanding  the  matrices  in  equation  (7.62),  we  get, 

x (tf )  =  $nx(t)  +  3?i2A(t)  +  Si  (*) 
A  (tf)  =  $2ix(t)  +  $22A(t)  +  g2  (t) 


(7.69) 


From  equation  (7.68), 

A  (tf)  =  TjiZ  (tf)  +  772X  (tf )  (7.70) 

Substituting  for  x(tf)  (in  the  above  equation)  from  Equation  (7.69),  we  get, 

A(tf)  =  ?ji$ux(t)  +  771^12 A (t)  +  7?xgi  (t)  +  772Z(tf)  (7.71) 

Substituting  the  above  equation  back  into  Equation  (7.69),  we  get, 

17l$lix(t)  +  7?l$12A(t)  +  77!gi(t)  +  772 Z  (tf)  =  $2ix(t)  +  $22A(t)  +  g2  (t) 

(7.72) 

By  rearranging  the  terms  in  the  above  equation,  we  get, 

[$22  -  771^12]  A  (t)  =  [771^11  -  $2i]x(t)  +  [77igi  (t)  -  g2  (t)  +  772Z  (tf)]  (7.73) 


k 
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From  equation  (7.73),  the  A(t)  needed  for  input  design  in  Equation  (7.56)  may 
be  derived  as, 

'  A(t)  =  [^22  —  X  (t)  +  [7j!gX  (t)  —  g3  (t)  +  (tf)]} 

(7.74) 

a  or, 

A(t)  =  L(t)x(t)  -  g(t)  (7.75) 

where  Z(t)  is  the  desired  output,  and  L  and  g  are  solutions  of  the  differential 
Riccatti  equations  shown  below: 

-  LHi  -  H?L  +  LH2L  4-  H4  =  L 

-  [LH2  -  Hf  ]  g  -  v2  +  LVx  =  g  (7.76) 

The  boundary  conditions  for  the  above  system  of  equations  are  given  as,  L(tf)  = 
and  g(tf)  =  7/2Z(tf),  which  are  constant  known  values.  With  these  boundary 
conditions,  solving  Equation  (7.76)  iteratively,  the  control  law  for  designing  the 
process  parameters  may  be  obtained. 

To  a  large  extent,  the  design  obtained  using  this  control  law  depends  on  the 
proper  selection  of  the  weighting  matrices  (Eq.  (7.48)).  Since  this  is  a  vital 
issue  in  getting  an  acceptable  design,  the  weighting  matrices  have  to  be  carefully 
designed  so  as  to  obtain  meaningful  inputs.  In  this  work,  the  weighting  matrices 
Q,R,  and  S  are  determined  as  follows  [69]: 

S-1  =  n  x  max.  of  diagonal  of  e  (tf)  eT  (tf), 

Q-1  =  n(tf  —  tQ)  x  max.  of  diagonal  of  e  (t)  eT  (t),  (7.77) 

4  m 

R-1  =  n/(tf  —  t0)  x  max.  of  diagonal  of  u  (t)  uT  (t). 
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where  e(tf)  is  known  from  the  desired  design  values.  The  values  of  e(t)  and  u(t) 
are  calculated  at  every  simulation  time  step  along  with  the  iterative  solution 
process. 

7.6  Performance  evaluation  and  selection  of  reduced  order  models 

The  order  of  the  full-size  state  space  model  was  reduced  using  the  method- 
s  discussed  in  earlier  sections  of  this  chapter.  A  comparative  study  was  then 
carried  out  to  evaluate  which  model  reduction  method  was  best  suited  for  metal 
forming  applications.  Two  numerical  examples  were  used  to  perform  the  compar¬ 
ative  study  of  the  model  reduction  methods  studied.  The  selection  of  a  suitable 
reduced  order  model  was  based  on  the  designed  die  velocities.  Of  the  two  exam¬ 
ples  presented  here,  the  first  example  has  no  friction  specified,  while  the  second 
example  uses  a  shear  friction  factor  of  0.3. 


7.6.1  H-block  compression  -  Frictionless  case 

A  closed  die  forging  was  simulated  to  deform  a  cylindrical  billet  into  a  H-block 
using  two  dies.  The  die  had  a  height  to  web  ratio  of  0.5.  A  cylindrical  billet  of 
radius  of  60  mm  and  height  40  mm  was  used  as  the  starting  workpiece.  Due  to  its 
symmetry,  only  one  half  of  the  workpiece  was  modeled.  The  discretized  workpiece 
had  96  elements  and  117  nodes  (Fig.  24).  The  simulation  was  carried  out  with 
an  initial  die  velocity  of  0.1  mm/sec  for  the  upper  die,  while  the  bottom  die  was 
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Figure  24:  FEM  discretization 
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kept  stationary.  A  rate  sensitive  billet  material  was  used  for  the  simulation.  The 
constitutive  (flow)  equation  for  the  material  is  given  as: 

a  =  k  flh  (7.78) 

where  k  is  a  proportionality  constant  used  in  the  flow  formulation,  g  is  the  rate 
sensitivity  factor,  and  h  is  the  effective  strain  exponent.  For  this  example  the 
values  of  k,  g,  and  h  were  10,  0.1,  and  0.0,  respectively.  The  strain  rate  of  a 
single  element  was  controlled,  and  element  number  53  was  chosen  as  the  element 
of  interest.  This  element  was  considered  critical  because  it  goes  through  a  number 
of  modes  of  deformation  during  the  simulation  due  to  its  location  in  the  billet. 
In  this  example  the  friction  factor  at  the  die/billet  interface  was  assumed  to  be 
zero. 

During  the  start  of  the  simulation,  seven  nodes  were  touching  the  (moving) 
die.  For  each  of  these  nodes,  the  normal  component  of  velocity  is  a  known 
quantity  and  is  not  considered  as  a  degree  of  freedom.  In  addition,  there  are 
four  nodes  associated  with  the  element  of  interest,  each  having  two  degrees  of 
freedom.  Since  the  degrees  of  freedom  of  the  system  are  also  considered  as  the 
states  during  state  space  modeling,  this  results  in  a  full-size  state  space  model 
with  15  states.  The  initial  strain  rate  of  element  53  was  0.0044/ sec  and  the 
desired  strain  rate  was  0.005/sec.  The  15  initial  states  were  reduced  to  4  using  the 
model  reduction  methods  described  earlier.  The  number  of  states  in  the  reduced 
model  was  decided  from  the  controllability  grammians  of  the  full  matrix.  In 
this  case,  only  four  grammians  had  significant  magnitude  and  the  rest  were  very 
small,  resulting  in  four  states  in  the  reduced  order  model.  Since  a  transformation 
matrix  was  used  to  construct  the  reduced  models,  the  weighting  matrices  were 
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modified  using  these  tansformation  matrices,  so  as  to  have  a  uniform  testing 
condition  for  all  the  models.  The  control  scheme  was  applied  for  designing  the 
die  velocity  using  the  state  space  models  derived  using  the  Aggregation,  Davison, 
Marshall,  Nicholson  and  BMR  approaches.  The  design  of  velocities  was  carried 
out  for  20  simulation  time  steps. 

The  comparative  results  are  shown  in  Fig.  25.  In  this  figure,  the  optimal  die 
velocities  are  plotted  as  a  function  of  time  for  20  time  steps  (0.4  sec).  The  simu¬ 
lations  were  conducted  using  a  time  step  of  0.02  sec.  At  each  time  step  the  state 
space  model  (A,  B  and  C)  is  built  based  on  the  current  billet  geometry,  die 
velocity  (designed  in  the  previous  step),  and  the  prevailing  frictional  conditions. 
For  the  first  step  the  die  velocity  is  taken  from  the  input  information  specified. 
A  finite  time  control  law  based  on  the  linear  quadratic  regulator  (LQR)  theory  is 
used  to  design  the  die  velocities  using  the  full  state  model.  These  results  are  used 
as  a  reference  against  which  the  results  of  the  reduced  order  model  are  compared. 
The  designed  velocities  (Fig.  25)  show  a  sharp  increase  in  the  initial  phase  of 
deformation  and  then  gradually  stabilize  with  time  as  the  strain  rate  approach¬ 
es  the  desired  value.  From  the  full  model,  the  earlier  mentioned  reduced-order 
systems  are  generated  and  an  optimal  die  velocity  is  again  designed.  It  may  be 
observed  that  the  velocities  designed  using  the  BMR  model  match  very  closely 
those  of  the  full  model.  The  maximum  percentage  error  in  the  design  was  0.3%. 
This  percentage  error  was  observed  in  the  initial  stages  of  the  velocity  design, 
and  the  error  decreased  as  the  simulation  progressed.  In  the  Aggregation  model, 
the  initial  results  (first  three  steps)  followed  the  trend  of  the  full  model.  But, 
after  that  the  designed  velocities  displayed  a  trend  quite  different  from  that  of 
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Figure  25:  Comparison  of  velocities  using  reduced  models  and 
full  model 
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the  full  model.  This  was  probably  because  the  transformation  (during  reduction) 
did  not  retain  the  controllable  characteristics  of  the  full  model.  Further,  in  the 
case  of  the  Aggregation  model,  the  performance  of  the  reduced  model  improved 
as  the  number  of  retained  states  was  increased.  Davison’s  model  resulted  m  a 
velocity  trend  similar  to  that  of  the  full  model  except  at  the  last  time  step  of 
the  simulation.  In  this  case  the  designs  do  not  match  because  only  the  first  four 
modes  are  retained  in  the  reduced  system.  Though  the  full  model  is  compensated 
for  within  the  reduced  system,  the  initial  response  of  the  unretained  states  may 
result  in  a  mismatch  in  the  results.  The  results  of  Marshall’s  and  Nicholsons’s 
method  are  similar,  and  to  a  large  extent  the  trend  of  the  full  state  model  is 
maintained,  in  both  cases.  Because  these  models  are  derived  from  model  trans¬ 
formations,  the  numerical  correspondence  with  the  full  model  is  not  exact.  In 
these  methods  the  proper  selection  and  grouping  of  the  dominant  modes  of  the 
system  largely  influence  the  system  characteristics  and  output. 

For  the  first  example,  the  reduced  model  built  using  the  BMR  scheme  retains 
the  controllable  and  observable  characteristics  of  the  full  model.  As  a  result,  the 
velocities  designed  using  the  BMR  model  and  full  model  match  very  closely. 
The  velocities  designed  using  the  Aggregation,  BMR,  Marshall,  Davison  and 
Nicholson  methods  along  with  results  of  the  full  model  are  tabulated  in  Tables 
(3)  through  (7),  respectively.  The  forging  simulation  was  earned  out  using  these 
velocities  and  the  finite  element  model  described  earlier.  The  resulting  strain 
rate  trajectories  obtained  are  depicted  in  Fig.  (26). 


jt 


135 


Table  4.  Comparison  of  optimized  velocities  for  BMR  model  (Frictionless  case) 


p 


Step 

Time 

Full  Model 

Red.  Model 

% 

No. 

(sec) 

(mm/sec) 

(mm/sec) 

error 

1 

0.02 

-0.2497 

-0.2504 

-0.2803 

2 

0.04 

-0.2387 

-0.2393 

-0.2514 

3 

0.06 

-0.2333 

-0.2338 

-0.2143 

4 

0.08 

-0.2288 

-0.2292 

-0.1748 

5 

0.10 

-0.2252 

-0.2255 

-0.1332 

6 

0.12 

-0.2223 

-0.2225 

-0.0900 

7 

0.14 

-0.2200 

-0.2202 

-0.0909 

8 

0.16 

-0.2182 

-0.2184 

-0.0917 

9 

0.18 

-0.2169 

-0.2170 

-0.0461 

10 

0.20 

-0.2160 

-0.2160 

0.0000 

11 

0.22 

-0.2153 

-0.2153 

0.0000 

12 

0.24 

-0.2150 

-0.2149 

0.0465 

13 

0.26 

-0.2149 

-0.2148 

0.0465 

14 

0.28 

-0.2150 

-0.2149 

0.0465 

15 

0.30 

-0.2153 

-0.2151 

0.0929 

16 

0.32 

-0.2157 

-0.2155 

0.0927 

17 

0.34 

-0.2163 

-0.2161 

0.0925 

18 

0.36 

-0.2169 

-0.2167 

0.0922 

19 

0.38 

-0.2177 

-0.2175 

0.0919 

20 

0.40 

-0.2186 

-0.2183 

0.1372 
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Table  5 


Comparison  of  optimized  velocities  for  Marshall’s  method  (Frictionless  case) 


% 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.1069 

0.1787 

-67.1656 

2 

0.04 

0.1102 

0.1846 

-67.5136 

3 

0.06 

0.1137 

0.1908 

-67.8100 

4 

0.08 

0.1173 

0.1973 

-68.2012 

5 

0.10 

0.1211 

0.2041 

-68.5384 

6 

0.12 

0.1251 

0.2111 

-68.7450 

7 

0.14 

0.1293 

0.2186 

-69.0642 

8 

0.16 

0.1337 

0.2263 

-69.2595 

9 

0.18 

0.1384 

0.2345 

-69.4364 

10 

0.20 

0.1433 

0.2431 

-69.6441 

11 

0.22 

0.1485 

0.2521 

-69.7643 

12 

0.24 

0.1540 

0.2616 

-69.8701 

13 

0.26 

0.1599 

0.2715 

-69.7936 

14 

0.28 

0.1661 

0.2821 

-69.8374 

15 

0.30 

0.1727 

0.2932 

-69.7742 

16 

0.32 

0.1797 

0.3049 

-69.6717 

17 

0.34 

0.1873 

0.3173 

-69.4074 

18 

0.36 

0.1953 

0.3305 

-69.2268 

19 

0.38 

0.2039 

0.3445 

-68.9554 

20 

0.40 

0.2132 

0.3594 

-68.5741 
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Table  6 


Comparison  of  optimized  velocities  for  Davison’s  method  (Frictionless  case) 


» 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.1069 

0.1318 

-23.2928 

2 

0.04 

0.1102 

0.1408 

-27.7677 

3 

0.06 

0.1137 

0.1504 

-32.2779 

4 

0.08 

0.1173 

0.1605 

-36.8287 

5 

0.10 

0.1211 

0.1712 

-41.3708 

6 

0.12 

0.1251 

0.1825 

-45.8833 

7 

0.14 

0.1293 

0.1945 

-50.4254 

8 

0.16 

0.1337 

0.2072 

-54.9738 

9 

0.18 

0.1384 

0.2207 

-59.4653 

10 

0.20 

0.1433 

0.2350 

-63.9916 

11 

0.22 

0.1485 

0.2503 

-68.5522 

12 

0.24 

0.1540 

0.2666 

-73.1169 

13 

0.26 

0.1599 

0.2841 

-77.6735 

14 

0.28 

0.1661 

0.3028 

-82.2998 

15 

0.30 

0.1727 

0.3229 

-86.9716 

16 

0.32 

0.1797 

0.3445 

-91.7084 

17 

0.34 

0.1873 

0.3679 

-96.4229 

18 

0.36 

0.1953 

0.3933 

-101.3825 

19 

0.38 

0.2039 

0.4210 

-106.4738 

20 

0.40 

0.2132 

0.4512 

-111.6323 
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Table  7 


Comparison  of  optimized  velocities  for  Nicholson’s  method  (Frictionless  case) 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.1069 

0.1318 

-23.2928 

2 

0.04 

0.1102 

0.1408 

-27.7677 

3 

0.06 

0.1137 

0.1504 

-32.2779 

4 

0.08 

0.1173 

0.1605 

-36.8287 

5 

0.10 

0.1211 

0.1712 

-41.3708 

6 

0.12 

0.1251 

0.1825 

-45.8833 

7 

0.14 

0.1293 

0.1945 

-50.4254 

8 

0.16 

0.1337 

0.2072 

-54.9738 

9 

0.18 

0.1384 

0.2207 

-59.4653 

10 

0.20 

0.1433 

0.2350 

-63.9916 

11 

0.22 

0.1485 

0.2503 

-68.5522 

12 

0.24 

0.1540 

0.2666 

-73.1169 

13 

0.26 

0.1599 

0.2841 

-77.6735 

14 

0.28 

0.1661 

0.3028 

-82.2998 

15 

0.30 

0.1727 

0.3229 

-86.9716 

16 

0.32 

0.1797 

0.3445 

-91.7084 

17 

0.34 

0.1873 

0.3679 

-96.4229 

18 

0.36 

0.1953 

0.3933 

-101.3825 

19 

0.38 

0.2039 

0.4210 

-106.4738 

20 

0.40 

0.2132 

0.4512 

-111.6323 
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7.6.2  H-block  compression  -  Friction  case 


0 


In  this  example,  a  closed  die  forging  of  an  H-block  from  a  cylindrical  billet  was 
simulated  using  one  die.  The  billet  considered  for  simulation  had  a  height  of  2.45 
inches  and  a  radius  of  2.45  inches.  Due  to  the  existing  symmetric  conditions,  only 
a  quarter  of  the  die  and  workpiece  were  modeled.  During  modeling  the  workpiece 
was  discretized  into  240  nodes  and  203  elements  (Fig.  27).  The  die  had  a  height  to 
web  ratio  of  1.0.  The  same  material  as  used  in  the  previous  example  was  used  for 
this  simulation.  The  frictional  condition  at  the  billet-die  interface  was  specified  by 
using  a  constant  shear  friction  factor  of  0.3.  The  simulation  was  initiated  with  an 
initial  die  velocity  of  0.1  in/sec.  The  control  of  a  single  element  was  considered  and 
a  desired  strain  rate  of  0.05 j  sec  in  element  number  153  was  taken  as  the  control 
objective.  This  element  is  critical  because  it  is  adjacent  to  the  die-billet  contact 
surface,  and  is  sensitive  to  the  frictional  conditions  at  the  die-billet  interface.  During 
the  first  step  of  simulation,  25  nodes  were  touching  the  die.  Each  of  these  nodes 
has  one  degree  of  freedom  as  explained  in  the  previous  example.  In  addition,  the 
element  of  interest  is  associated  with  four  nodes  having  two  degrees  of  freedom 
each.  This  results  in  a  full-size  (condensed)  state  space  system  with  33  states.  At 
the  end  of  the  second  time  step,  29  nodes  were  touching  the  die.  As  explained 
above,  this  would  result  in  a  state  space  model  having  37  states.  By  applying  the 
model  reduction  techniques  explained  earlier,  the  number  of  states  was  reduced 
to  10.  In  this  case,  the  number  of  states  in  the  reduced  model  was  decided  by 
the  controllability  grammians.  In  the  current  example  the  first  10  grammians  had 
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significant  values  compared  to  the  rest,  resulting  in  a  reduced  order  model  having 
10  states. 

To  evaluate  the  perormance  of  the  reduced  order  models,  the  input  velocities 
were  designed  using  the  full  model  and  compared  with  those  of  the  reduced  models. 
Because  the  strain  rate  desired  was  very  low  compared  to  the  initial  strain  rate  of 
element  153,  the  design  trend  using  the  full  model  revealed  that  the  die  velocities 
decreased  rapidly  as  the  simulation  progressed.  The  consolidated  results  and  com¬ 
parisons  of  the  designed  velocities  are  depicted  in  Fig.  (28).  The  simulations  were 
carried  out  using  a  time  step  of  0.02  sec  for  20  time  steps.  It  may  be  observed  that 
the  results  of  the  BMR  model  match  almost  exactly  with  that  of  the  full  model 
(Table  8),  with  a  maximum  error  of  0.06%.  The  designed  ram  velocities  using  the 
Aggregation,  Marshall,  Davison  and  Nicholson  methods  were  also  compared  with 
those  of  the  full  model.  These  results  are  tabulated  in  Tables  (9)  through  (12), 
respectively.  Here  again,  the  Aggregation  model  performed  well  in  the  initial  few 
steps.  But,  after  the  seventh  step  the  results  from  this  model  did  not  conform  with 
that  of  the  full  model.  This  behavior  of  the  aggregated  model  is  due  to  the  trun¬ 
cation  of  the  higher  eigenvalues  of  the  system  during  the  reduction  process.  The 
other  models  showed  trends  in  close  conformance  with  the  full-size  model.  It  is 
very  likely  that  retaining  more  states  in  the  reduced  models  would  result  in  more 
accurate  results. 


* 
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Table  8.  Comparison  of  optimized  velocities  for  BMR  model  (Friction  case) 


* 


Step 

No. 

Time 

(sec) 

Full  Model 
(in/sec) 

Red.  Model 
(in/sec) 

% 

error 

1 

0.02 

-0.002341 

-0.002342 

0.0548 

2 

0.04 

-0.009623 

-0.009629 

0.0570 

3 

0.06 

-0.016222 

-0.016229 

0.0489 

4 

0.08 

-0.022782 

-0.022792 

0.0423 

5 

0.10 

-0.029304 

-0.029315 

0.0365 

6 

0.12 

-0.035785 

-0.035796 

0.0313 

7 

0.14 

-0.042225 

-0.042236 

0.0266 

8 

0.16 

-0.048623 

-0.048634 

0.0223 

9 

0.18 

-0.054980 

-0.054990 

0.0184 

10 

0.20 

-0.061294 

-0.061304 

0.0149 

11 

0.22 

-0.067566 

-0.067574 

0.0117 

12 

0.24 

-0.073794 

-0.073800 

0.0088 

13 

0.26 

-0.079979 

-0.079984 

0.0061 

14 

0.28 

-0.086119 

-0.086122 

0.0038 

15 

0.30 

-0.092215 

-0.092217 

0.0016 

16 

0.32 

-0.098266 

-0.098266 

0.0003 

17 

0.34 

-0.104272 

-0.104270 

-0.0020 

18 

0.36 

-0.110232 

-0.110229 

-0.0035 

19 

0.38 

-0.116147 

-0.116141 

-0.0048 

20 

0.40 

-0.122016 

-0.122008 

-0.0060 
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Table  9 


Comparison  of  optimized  velocities  for  aggregation  method  (Friction  case) 


% 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.0993 

0.2760 

-177.9456 

2 

0.04 

0.1037 

0.2764 

-166.5381 

3 

0.06 

0.1083 

0.2767 

-155.4940 

4 

0.08 

0.1133 

0.2771 

-144.5719 

5 

0.10 

0.1186 

0.2775 

-133.9798 

6 

0.12 

0.1243 

0.2779 

-123.5720 

7 

0.14 

0.1303 

0.2782 

-113.5073 

8 

0.16 

0.1368 

0.2786 

-103.6550 

9 

0.18 

0.1436 

0.2790 

-94.2897 

10 

0.20 

0.1510 

0.2794 

-85.0331 

11 

0.22 

0.1588 

0.2797 

-76.1335 

12 

0.24 

0.1671 

0.2801 

-67.6242 

13 

0.26 

0.1760 

0.2805 

-59.3750 

14 

0.28 

0.1855 

0.2809 

-51.4286 

15 

0.30 

0.1956 

0.2813 

-43.8139 

16 

0.32 

0.2064 

0.2817 

-36.4825 

17 

0.34 

0.2179 

0.2821 

-29.4631 

18 

0.36 

0.2302 

0.2824 

-22.6759 

19 

0.38 

0.2434 

0.2828 

-16.1873 

20 

0.40 

0.2574 

0.2832 

-10.0233 

* 
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Table  10.  Comparison  of  optimized  velocities  for  Marshall’s  method  (Friction 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.0993 

0.2039 

-105.3374 

2 

0.04 

0.1037 

0.2045 

-97.2035 

3 

0.06 

0.1083 

0.2050 

-89.2890 

4 

0.08 

0.1133 

0.2055 

-81.3769 

5 

0.10 

0.1186 

0.2060 

-73.6931 

6 

0.12 

0.1243 

0.2066 

-66.2108 

7 

0.14 

0.1303 

0.2071 

-58.9409 

8 

0.16 

0.1368 

0.2076 

-51.7544 

9 

0.18 

0.1436 

0.2081 

-44.9164 

10 

0.20 

0.1510 

0.2087 

-38.2119 

11 

0.22 

0.1588 

0.2092 

-31.7380 

12 

0.24 

0.1671 

0.2097 

-25.4937 

13 

0.26 

0.1760 

0.2103 

-19.4886 

14 

0.28 

0.1855 

0.2108 

-13.6388 

15 

0.30 

0.1956 

0.2113 

-8.0266 

16 

0.32 

0.2064 

0.2119 

-2.6647 

17 

0.34 

0.2179 

0.2124 

2.5241 

18 

0.36 

0.2302 

0.2129 

7.5152 

19 

0.38 

0.2434 

0.2135 

12.2843 

20 

0.40 

0.2574 

0.2140 

16.8609 

Table  11.  Comparison  of  optimized  velocities  for  Davison’s  method  (Friction  case) 


Full  model 
(mm/sec) 
0.0993 
0.1037 
0.1083 
0.1133 
0.1186 
0.1243 
0.1303 
0.1368 
0.1436 
0.1510 
0.1588 
0.1671 
0.1760 
0.1855 
0.1956 
0.2064 
0.2179 
0.2302 
0.2434 
0.2574 


Red.  Model 
(mm/sec) 
0.1233 
0.1234 
0.1235 
0.1236 
0.1237 
0.1238 
0.1240 
0.1241 
0.1242 
0.1243 
0.1244 
0.1245 
0.1247 
0.1248 
0.1249 
0.1250 
0.1251 
0.1253 
0.1254 
0.1255 


-24.1692 

-18.9971 

-14.0351 

-9.0909 

-4.3002 

0.4023 

4.8350 

9.2836 

13.5097 

17.6821 

21.6625 

25.4937 

29.1477 

32.7224 

36.1452 

39.4380 

42.5883 

45.5691 

48.4799 

51.2432 


Table  12 


Comparison  of  optimized  velocities  for  Nicholson’s  method  (Friction  case) 


¥ 


Step 

No. 

Time 

(sec) 

Full  model 
(mm/sec) 

Red.  Model 
(mm/sec) 

% 

Error 

1 

0.02 

0.0993 

0.1233 

-24.1692 

2 

0.04 

0.1037 

0.1234 

-18.9971 

3 

0.06 

0.1083 

0.1235 

-14.0351 

4 

0.08 

0.1133 

0.1236 

-9.0909 

5 

0.10 

0.1186 

0.1237 

-4.3002 

6 

0.12 

0.1243 

0.1238 

0.4023 

7 

0.14 

0.1303 

0.1240 

4.8350 

8 

0.16 

0.1368 

0.1241 

9.2836 

9 

0.18 

0.1436 

0.1242 

13.5097 

10 

0.20 

0.1510 

0.1243 

17.6821 

11 

0.22 

0.1588 

0.1244 

21.6625 

12 

0.24 

0.1671 

0.1245 

25.4937 

13 

0.26 

0.1760 

0.1247 

29.1477 

14 

0.28 

0.1855 

0.1248 

32.7224 

15 

0.30 

0.1956 

0.1249 

36.1452 

16 

0.32 

0.2064 

0.1250 

39.4380 

17 

0.34 

0.2179 

0.1251 

42.5883 

18 

0.36 

0.2302 

0.1253 

45.5691 

19 

0.38 

0.2434 

0.1254 

48.4799 

20 

0.40 

0.2574 

0.1255 

51.2432 

* 
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The  reduced  model  constructed  using  the  BMR  method  resulted  in  a  velocity 
schedule  very  close  to  that  of  the  full  state  model  (Table  8),  and  the  corresponding 
strain  rate  profile  is  depicted  in  Fig.  29. 

7.7  Analysis  of  results 

In  this  work,  an  efficient  model  reduction  scheme  suitable  for  the  metal  forming 
applications  is  identified.  The  methodology  developed  is  general  purpose  and  ap¬ 
plicable  to  most  metal  forming  processes  such  as  forging,  extrusion,  upsetting  and 
coining  processes.  Of  the  possible  order  reduction  methods  studied,  the  balanced 
model  reduction  (BMR)  method  was  the  most  accurate  scheme  suitable  for  process 
control  of  metal  forming  operations.  Because  the  BMR  scheme  is  based  on  the  con¬ 
trollable  and  observable  properties  of  the  original  system,  the  designed  velocities 
matched  almost  exactly  with  those  of  the  full  model.  Another  advantage  of  this 
method  is  that  it  can  be  used  for  a  very  large  reductions  in  the  number  of  states. 
Also,  this  method  does  not  need  any  additional  calculation  of  weighting  matrices 
since  it  is  based  on  the  controllability  of  the  original  system.  This  results  in  a 
savings  in  computational  time  during  simulation,  as  compared  to  the  other  model 
reduction  methods. 

The  consolidated  results  are  presented  in  Figures  (25)  and  (29).  In  the  first  case, 
the  discretized  system  resulted  in  234  states,  which  was  reduced  to  15  states  by  using 
finite  element  condensation  techniques.  This  was  further  reduced  to  four  states  by 
the  model  reduction  schemes  considered.  The  initial  strain  rate  of  the  uncontrolled 
system  was  0.0044 /sec  and  the  objective  was  to  increase  this  to  0.005/sec.  Because 
the  desired  strain  rate  is  higher  than  the  initial  strain  rate,  the  control  dictated  a 
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Figure  29:  Comparison  of  desired  and  controlled  strain  rates 
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sharp  increase  in  the  velocity.  The  change  in  the  velocity  gradually  reduced  with  die 
stroke,  maintaining  the  strain  rate  at  0.005/sec.  This  trend  was  exhibited  by  all  the 
reduced  models  except  the  Aggregated  model.  The  reduced  model  constructed  using 
the  balanced  model  reduction  scheme  resulted  in  a  design  closest  (comparatively) 
to  that  of  the  the  full  model.  For  the  BMR  model,  the  maximum  deviation  in 
velocities  from  the  original  (full-size  model)  design  was  0.3%. 

In  the  second  case,  the  finite  element  model  had  240  nodes,  and  the  number  of 
states  in  the  original  system  was  480.  The  condensation  scheme  reduced  this  large 
number  of  states  to  37  at  the  second  step,  and  model  reduction  further  reduced  the 
number  of  states  to  10.  The  initial  strain  rate  of  element  153  was  0.0996/sec,  and 
the  desired  effective  strain  rate  was  0.005/sec.  Because  a  reduction  in  effective  strain 
rate  is  desired,  the  designed  velocities  decreased  sharply  to  ensure  that  the  desired 
strain  rate  is  achieved.  This  tendency  was  shown  by  all  the  reduced  models  except 
the  Aggregation  model.  A  characteristic  feature  of  the  modal  analysis  methods  is 
that  the  steady  state  response  matches  with  the  parent  model.  Hence,  the  reduced 
model  behaves  similar  to  the  full  model  after  the  settling  time  is  reached.  But  in 
metal  forming  processes,  the  simulation  time  step  cannot  be  large  because  of  the 
nonlinearity  of  the  system/process.  Therefore,  in  general,  modal  analysis  methods 
are  not  suitable  for  metal  forming  processes. 


CHAPTER  8 


Numerical  Examples  -  Reduced  Order  Models 

The  condensed  state  space  model  representing  the  coupled  thermomechamcal 
forging  system  is  first  constructed.  This  model  is  further  reduced  by  retaining  the 
controllable  and  observable  subspace  of  the  system  using  the  balanced  state  space 
representation.  The  process  parameters  are  then  designed  using  LQR  theory  as  an 
off-line  design  tool.  The  design  problem  is  posed  as  an  output  tracking  problem 
and  the  trajectory  to  be  tracked  is  the  desired  strain-rate  profile  of  the  element 
of  interest.  The  state  space  model  is  updated  at  the  end  of  each  time  step  to 
accomodate  the  frequent  changes  in  geometry  and  boundary  conditions. 

The  design  approach  is  demonstrated  using  two  example  cases.  In  the  first 
example,  the  design  of  process  parameters  (die  velocity)  is  performed  while  forging 
an  integrated  blade  and  rotor  (IBR)  disk  under  isothermal  conditions.  This  problem 
is  carried  out  with  two  different  initial  conditions.  In  the  first  case  an  initial  die 
velocity  of  1.0  in/sec  was  used,  and  in  the  second  case,  the  design  was  carried 
out  using  an  initial  ram  velocity  of  0.5  in/sec.  In  both  cases,  the  method  was 
demonstrated  using  two  different  strain  rate  requirements.  Also,  the  simulation 
was  carried  out  only  for  200  time  steps,  and  not  until  completion.  This  is  because 
the  goal  of  these  examples  is  limited  to  demonstrating  the  concept  developed  and 
not  for  simulating  and  analyzing  the  entire  forging  process.  In  the  second  example, 
die  velocity  design  was  performed  under  nonisothermal  conditions  for  an  engine  disk 


154 


forging.  For  validation  purposes  the  design  obtained  using  the  reduced  model  was 
compared  with  that  obtained  using  the  full  model. 


8.1  Isothermal  design  -  IBR  disk  forging 

Forging  of  an  IBR  disk  is  generally  conducted  in  multiple  stages  because  of  the 
complexity  and  intricacy  of  the  finished  part.  The  die/billet  shapes  for  this  forging 
were  determined  by  earlier  researchers  using  a  trial  and  error  approach.  Figure  (30) 
shows  the  billet  and  die  shapes  used  in  simulating  the  IBR  disk  forging.  Figures 
(30a)  and  (30b)  show  the  start  and  finish  of  the  first  stage  of  forging.  The  second 
stage  forging  simulations  are  shown  in  Figures  (30c)  and  (30d). 

In  this  work,  the  billet  and  die  shapes  shown  in  Figure  (31)  were  used  for  sim¬ 
ulation  while  designing  the  process  parameters.  Though  there  is  a  considerable 
change  in  billet  shape  during  the  first  stage  of  simulation,  this  process  is  similar 
to  pure  compression  where  only  one  mode  of  deformation  is  dominant.  On  the 
other  hand,  the  second  stage  has  multiple  modes  of  deformation  and  constitutes  a 
more  challenging  forging  problem.  As  a  result,  the  second  stage  of  the  IBR  disk 
forging  was  selected  for  demonstrating  the  methodology  proposed  in  this  work.  To 
evaluate  the  performance  of  the  reduced  order  model,  the  same  case  was  studied 
with  two  different  initial  ram  velocities  and  with  two  different  design  (strain  rate) 
requirements.  Element  275  was  chosen  as  the  element  of  interest  for  this  example 
(Figure  31).  This  element  is  ‘critical’  because  it  goes  through  a  variety  of  defor¬ 
mation  modes  as  the  forging  simulation  progresses.  A  hypothetical  rate  dependent 
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Figure  31:  Finite  eienaent  model  of  Integrated  Blade  and  Rotor 
disk 
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material  with  constitutive  equation  cr  =  ke9eh  was  used  in  the  simulation  [Ref.  Eqn. 
(7.78)].  In  this  case,  the  proportionality  constant  k ,  the  strain  rate  exponent  g ,  and 
the  strain  sensitivity  factor  h,  were  assumed  to  be  10.0,  0.1,  and  0,  respectively. 
Due  to  the  symmetry  of  the  model  only  one  half  of  the  disk  was  modeled.  The 
discretized  workpiece  had  302  quadrilateral  elements  and  347  nodes.  A  constant 
temperature  of  1700°F  was  maintained  in  the  workpiece  for  the  entire  isothermal 
forging  process.  The  bottom  die  was  kept  stationary  while  the  top  die  was  allowed 
to  move.  Interface  frictional  conditions  were  enforced  using  a  constant  shear  friction 
factor  of  0.15. 

In  the  first  case  the  simulation  was  started  with  an  initial  die  velocity  of  1.0 
in/sec.  Two  different  strain-rate  requirements  of  0.2/sec  and  0.5/sec  were  imposed 
on  element  275.  The  0.5/sec  strain  rate  case  was  compared  with  the  full  state  space 
model  results  and  the  designed  velocities  for  200  simulation  time  steps  are  shown 
in  Figure  (32).  The  full  model  results  were  used  as  the  reference  design  values. 
From  the  figure  it  can  be  observed  that  the  full  model  and  reduced  model  results 
match  exactly,  and  the  difference  in  results  obtained  from  the  two  models  is  zero. 
At  the  beginning  of  the  process,  the  full  state  space  model  had  26  states  and  the 
controllable  subspace  (reduced  order  model)  had  4  states.  At  the  completion  of  200 
time  steps,  there  were  34  states  in  the  full  model  and  5  states  in  the  reduced  model. 
Figure  (33)  depicts  the  effective  strain  rate  profile  of  element  275  for  the  the  full  and 
reduced  models.  The  difference  in  results  for  the  two  models  is  again  zero,  because 
the  balanced  model  consists  of  the  controllable  subspace  of  the  full  system  and  the 
weighting  matrices  (for  control)  used  in  both  the  simulations  are  the  same.  The 
fluctuation  in  die  velocity /strain  rate  towards  the  end  of  the  200i/l  step  is  because 
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of  the  change  in  boundary  conditions  of  the  system  as  nodes  get  detached  and 
reattached  to  the  die  surface.  With  the  same  initial  die  velocity,  a  different  design 
requirement  of  0.2/sec  was  now  imposed  on  element  275.  The  design  procedure  was 
then  repeated  for  200  simulation  steps.  The  strain  rate  of  element  275  obtained 
using  the  designed  ram  velocity  schedule  was  compared  with  that  obtained  using  a 
constant  velocity  of  1.0  in/sec  in  Figure  (34).  This  figure  shows  a  plot  of  the  strain 
rate  against  die  stroke  upto  a  stroke  length  of  0.0993  in. 

The  design  was  continued  for  the  strain-rate  requirement  of  0.5/sec,  until  ele¬ 
ments  distorted  significantly  and  remesh  was  needed.  The  reduced  model  at  this 
stroke  (0.98  in.)  had  5  states,  while  the  full  model  had  66  states.  The  velocity 
and  strain  rate  schedules  for  this  example  are  described  in  Figure  (35)  and  Figure 
(36),  respectively.  It  may  be  observed  that  there  is  a  sudden  jump  in  strain-rate 
(Figure  (37))  at  a  stroke  of  0.17  in.  due  to  the  detaching  of  a  node  from  the  bottom 
die.  Further,  a  series  of  oscillations  in  strain-rate  may  be  observed  at  a  stroke  of 
0.4  in.  These  oscillations  may  be  attributed  to  the  attaching  of  additional  nodes 
to  the  die  at  this  stroke.  As  new  nodes  come  in  contact  with  the  die,  the  time 
step  used  in  ALPID  decreases  considerably.  The  controller  now  has  a  smaller  time 
duration  to  meet  the  strain  rate  requirement  in  the  ‘control’  element.  It  thus  makes 
large  corrections  in  the  input  (die  velocity).  If  these  adjustments  are  not  within  the 
equipment  acceleration  limit,  the  die  velocity  tends  to  bounce  back  and  forth  be¬ 
tween  the  acceleration  bounds  resulting  in  the  fluctuations  observed  in  Figure  (37). 
For  on-line  control  implementation,  these  oscillations  must  be  smoothed  about  their 
mean  value  at  the  disturbance  locations  and  then  used. 
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Figure  34:  Comparison  of  strain  rates  for  IBR  forging  using  reduced 
order  model 
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Figure  35:  Die  velocity  profiles  -  Designed  and  constant  velocity 
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Figure  36:  Strain  rate  profiles  -  Designed  and  constant  velocity 
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Figure  37:  Comparison  of  designed  velocities  and  contact  nodes  vs. 
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Number  of  nodes 


In  the  second  case,  an  initial  die  velocity  of  0.5  in/sec  was  used  with  two  d- 
ifferent  design  requirements  of  0.15/sec  and  0.1/sec  on  the  strain  rate  of  element 
275.  Figure  (38)  shows  a  comparison  between  the  strain  rate  profiles  for  the  con¬ 
trolled  and  uncontrolled  (constant  velocity)  cases  for  200  simulation  steps.  The 
figure  demonstrates  the  effectiveness  of  the  reduced  order  model  in  meeting  the  de¬ 
sign  requirement  under  different  initial  conditions  and  confirms  the  stability  of  the 
proposed  approach. 


8.2  Non-isothermal  design  -  Engine  disk  forging 

In  the  second  example,  process  parameter  design  for  an  engine  disk  forging 
was  carried  out  under  nonisothermal  conditions.  The  finite  element  discretization 
of  the  billet  and  die  for  the  engine  disk  are  shown  in  Figure  (39).  An  initial  die 
temperature  of  600°F  and  initial  billet  temperature  of  1700°F  was  chosen  for  this 
simulation.  Figure  (39)  also  shows  the  location  of  the  element  of  interest  (element 
23)  and  critical  node  (node  18)  for  this  example.  The  critical  node  may  be  defined 
as  the  node  at  the  die-billet  interface  whose  temperature  is  to  be  tracked.  Only 
the  interface  nodes  qualify  to  be  critical  nodes  because  they  go  through  a  variety 
of  heat  transfer  modes  as  deformation  proceeds.  The  design  requirement  was  to 
maintain  the  strain-rate  of  the  element  23  at  0.7/sec.  Only  a  quarter  of  the  billet 
was  modeled  because  of  the  inherent  symmetry  of  the  engine  disk  along  the  two 
cartesian  axes.  The  billet  was  discretized  using  50  quadrilateral  elements  and  66 
nodes.  The  billet  material  used  for  this  example  was  titanium  alloy  (Ti  6242),  and 
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the  design  was  started  using  a  single  die  (top  die)  moving  with  a  velocity  of  1.0 
in/sec. 

The  full-size  state  space  model  for  this  system  (nonisothermal)  resulted  in  21 
states.  The  number  of  states  was  further  reduced  to  10  by  using  the  BMR  method. 
Process  parameter  (die  velocity)  design  was  performed,  and  Figure  (40)  shows  the 
resulting  strain  rate  trajectory.  In  the  figure,  the  results  from  the  full  state  model 
were  used  as  the  reference.  It  may  be  observed  that  the  full  and  reduced  order 
models  gave  exactly  the  same  results.  The  corresponding  die  velocity  profiles  for 
the  two  models  (Figure  (41))  also  matched  exactly.  The  temperature  profile  of 
node  18  was  compared  for  the  two  models  and  is  depicted  in  Figure  (42).  Again,  it 
was  observed  that  the  results  of  the  full-size  and  reduced  order  state  space  models 
match  very  well,  with  a  maximum  difference  in  nodal  temperature  of  0.25%.  The 
simulation  was  carried  out  till  the  die  was  completely  filled,  and  at  the  end  of  the 
simulation  it  was  observed  that  the  full  model  had  41  states  while  the  reduced  order 
model  had  a  total  of  11  states. 
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Figure  39:  Finite  element  model  of  die  and  billet  for  engine  disk 
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Figure  40:  Comparison  of  strain  rates  for  engine  disk  forging 
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Figure  42:  Comparison  of  nodal  temps,  (node  18)  for  engine  disk 
forging 
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CHAPTER  9 


Practical  Aspects  in  Process  Control  Design 


9.1  Smoothened  die  velocity  profiles 

The  optimal  control  scheme  described  earlier  worked  satisfactorily  for  most 
metal  forming  problems.  But,  it  was  observed  that  simulations  involving  intricate 
and  complex  die/billet  shapes  resulted  in  frequent  changes  in  the  material  flow 
direction  and  boundary  conditions  of  the  system.  This  implies  that  nodes  get 
attached/detached  from  the  die  boundary  at  frequent  intervals  as  the  simulation 
progresses.  This  phenomemon  causes  the  optimal  control  scheme  to  produce  an 
oscillating  die  velocity  schedule.  In  addition,  there  is  a  complex  interplay  of 
forces  during  die-fill  which  results  in  unsteady  and  non-uniform  material  flow. 
This  again  causes  the  control  algorithm  to  produce  a  fluctuating  die  velocity 
during  the  final  stages  of  the  forming  process.  A  fluctuating  or  oscillating  die 
velocity  schedule  is  impractical  to  implement  on  a  hydraulic  press.  Therefore, 
designed  velocity  curves  were  smoothened  locally  using  curve-fitting  techniques 
before  implementation.  To  check  the  effect  of  this  approximation  on  the  results 
obtained,  the  new  velocity  schedule  was  fed  back  into  the  system  using  the  same 
design  conditions  and  requirements  as  the  original  (control)  simulation.  Fig.  43 
shows  a  comparison  of  the  original  (designed)  and  smoothened  velocity  profiles 
for  an  axisymmetric  engine  disk  forging.  The  corresponding  strain  rate  profiles 
are  depicted  in  Fig.  44.  It  is  observed  that  the  strain  rate  obtained  using  the 
smoothened  velocity  profile  satisfies  the  (strain  rate)  design  requirement,  and  is 
also  in  close  conformance  with  that  obtained  using  the  original  velocity  profile. 


A  smoothened  velocity  schedule  was  also  produced  for  a  plane  strain  channel 
section  forging.  Again,  it  was  observed  that  the  strain  rate  obtained  met  the 
design  requirement  satisfactorily,  and  was  in  conformance  with  that  obtained 
using  the  original  velocity.  Figures  45  and  46  depict  the  die  velocity  and  strain 
rate  profiles  for  this  example  case.  It  may  thus  be  concluded  that  the  smoothened 
velocity  curves  are  an  acceptable  approximation,  and  may  be  used  for  practical 
application  of  the  control  algorithm  described  in  this  work. 

In  addition  to  the  advantage  during  implementation,  the  smoothened  velocity 
profiles  are  also  an  asset  while  carrying  out  parametric  studies  for  metal  forming 
simulations.  In  such  studies,  generally,  repetitive  simulations  must  be  performed, 
and  previously  constructed  smoothened  velocity  profiles  can  directly  be  used  for 
carrying  out  additional  simulations,  instead  of  using  the  computationally  inten¬ 
sive  control-based  design  methodology.  This  would  save  a  considerable  amount  of 
computational  time  and  effort  during  the  simulation,  since  expensive  numerical 
calculations  involved  in  the  design  process  may  be  avoided. 

9.2  Stability  and  validity  of  the  state  space  models 

The  results  obtained  using  the  control-based  methodology  depend,  to  a  large 
extent,  on  the  prevailing  properties  and  characteristics  of  the  metal  forming 
system.  This  means  that  the  designed  velocities  and  the  effectiveness  of  the 
control  scheme  are  dependent  on  the  system  matrices  (plant,  input,  and  output 
matrices)  and  the  weighting  matrices  used  in  the  LQR  design  scheme.  But,  for 
simple  forging  simulations,  having  uniform  and  well  distributed  metal  flow,  it  was 
observed  that  the  system  matrices  did  not  change  significantly  from  one  time  step 
to  the  next.  Studies  were  thus  carried  out  to  determine  the  effect  of  using  only 
one  set  of  system  matrices  for  the  entire  simulation,  instead  of  building  the  state 
space  model  at  every  time  step.  This  would  again  reduce  the  computational  cost 
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Figure  43:  Engine  disk  forging  -  Velocity  profiles 
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Figure  44:  Engine 
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Figure  45:  Channel  section  forging  -  Velocity  profiles 
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Figure  46:  Channel  section  forging  -  Strain  rate  profiles 
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and  time  significantly.  An  axisymmetric  engine  disk  forging  was  considered  as 
the  test  case,  and  the  design  of  velocities  was  carried  out  for  the  following  three 
situations: 

1.  State  space  model  built  at  every  time  step. 

2.  State  space  model  built  only  once  at  the  beginning  of  the  simulation. 

3.  State  space  model  built  whenever  boundary  condition  changes. 

In  Case  3,  the  boundary  condition  changed  12  times  during  the  entire  sim¬ 
ulation,  which  lasted  about  130  time  steps.  As  a  result,  the  state  space  model 
had  to  be  built  12  times  during  the  design  process.  On  an  average,  each  state 
space  model  was  valid  for  about  10  simulation  time  steps  for  this  example.  Case 
3  is  also  more  accurate  than  Case  2  because  it  accounts  for  the  change  in  the 
size  of  the  system  matrices  as  deformation  progresses.  In  other  words,  in  case 
3,  the  system  matrices  are  built  and  used  (for  designing  the  die  velocity)  when¬ 
ever  nodes  attach/detach  from  the  die  surface.  Figures  47  and  48  show  the  die 
velocity  and  strain  rate  profiles  obtained  for  the  above  three  cases.  It  may  be 
observed  from  Fig.  47  that  the  die  velocities  designed  for  the  above  three  cases 
match  closely.  Fig.  48  shows  that  the  strain  rate  requirement  is  met  effectively 
for  all  the  three  cases.  Since  the  results  from  case  2  match  so  closely  with  that 
of  the  other  two  cases,  it  may  be  concluded  that  for  this  case  the  one  state  space 
model  built  at  the  beginning  of  the  simulattion  was  valid  for  the  entire  duration 
of  the  simulation.  In  addition,  use  of  approximate  system  matrices  for  the  entire 
simulation  has  a  smoothening  effect  on  the  designed  velocities,  thereby  reducing 
the  fluctuations  in  the  designed  die  velocities.  This  in  itself  is  a  big  advantage, 
as  explained  in  the  previous  section.  It  may  thus  be  concluded  that  for  simple 
forging  problems,  approximate  system  matrices  may  be  used  to  design  die  veloc¬ 
ities  without  radically  altering  the  performance  of  the  system.  In  addition,  by 


intelligently  selecting  points  along  the  simulation  path  where  the  system  matri¬ 
ces  are  to  be  built,  a  large  amount  of  computational  time  may  be  saved  without 
sacrificing  the  accuracy  of  the  system. 

The  above  argument  may  be  extended  and  applied  to  the  design  of  weighting 
matrices  also.  In  the  original  design  scheme,  the  weighting  matrices  for  control 
are  determined  at  every  time  step  of  the  simulation.  For  systems  behaving  in  a 
fairly  linear  and  predictable  manner,  constant  weighting  matrices  may  be  used 
for  the  entire  simulation  without  affecting  the  results  significantly.  Figs.  49  and 
50  show  the  velocities  and  (corresponding)  strain  rates  obtained  using  constant 
and  varying  weighting  matrices  for  an  engine  disk  forging.  In  both  cases,  the 
conformance  with  the  original  results  is  very  good.  Therefore,  it  may  be  con¬ 
cluded  that  the  use  of  constant  weighting  matrices  in  the  control  scheme  is  again 
an  acceptable  approximation. 

All  the  above  mentioned  approximations  made  in  the  formulation  and  imple¬ 
mentation  procedures  result  in  a  large  savings  in  computational  time  and  effort 
during  the  simulation,  while  still  realizing  effective  control  and  design.  However, 
it  must  be  noted  that  these  practical  aspects  are  related  to  the  simulation  and 
are  not  directly  applicable  to  the  actual  forging  process  itself. 
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Figure  47:  Die  velocities  using  constant  and  varying  system  matri- 
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Figure  48:  Strain  rates  using  constant  and  varying  system  matrices 
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Figure  49:  Die  velocities  using  constant  and  varying  weighting 
matrices 
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Figure  50:  Strain  rates  using  constant  and  varying  weighting  ma¬ 
trices  t 
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9.3  Eigen  mode  analysis  for  process  characterization 


This  section  presents  a  methodology  for  identifying  the  tension  and  compres¬ 
sion  zones  in  the  wrokpiece  material,  based  on  the  study  of  eigenmodes  of  the 
deforming  system.  The  process  of  identification  of  the  tension  and  compression 
zones  is  then  related  to  the  selection  of  favorable  operating  zones/conditions  for 
the  process,  and  for  predicting  the  stability  of  the  deforming  material. 


Elemental  level  eigenmode  analysis 

The  derivation  of  the  finite  element  governing  equations  for  metal  forming 
analysis  have  been  described  earlier  in  chapter  2  of  this  report.  In  the  finite 
element  formulation,  the  objective  is  to  design  admissible  velocity  fields  by  min¬ 
imizing  the  potential  energy  of  the  system  (Eq.  (2.8))  while  also  maintaining 
the  incompressibility  condition  on  the  velocity  field.  This  is  done  by  introducing 
a  penalty  constant  for  the  volumetric  component  and  converting  the  variational 
problem  into  a  stationary  value  problem,  whose  first  order  variation  (Eq.  (2.10)) 
is  stated  as, 

Stt  =  f  dSedV  +  K  [  e'vSeydV  —  f  F{8uidS  =  0  (9*1) 

Jv  JV  JsF 

where  K  is  the  penalty  parameter,  and  is  generally  a  large  positive  quantity. 

On  the  right  hand  side  of  Eq.  (9.1),  the  first,  second,  and  third  terms  rep¬ 
resent  the  shear  component,  the  volumetric  component,  and  the  work  done  by 


external  tractional  forces,  respectively.  Because  of  the  penalty  parameter,  the 
volumetric  deformation  component  is  easily  separable  from  the  shear  deformation 
component.  In  a  similar  manner  the  stiffness  matrix  of  the  element  can  also  be 
separated  into  the  deformation  part  and  volumetric  part.  The  stiffness  equations 
in  metal  forming  analysis  are  nonlinear  in  nature,  and  the  velocity  field  solution 
is  generally  obtained  in  an  iterative  manner.  At  the  converged  solution  point, 
the  stiffness  matrix  of  any  element  may  be  represented  as, 

Kc  =  Kd  +  Kyo  (9-2) 

where  Kd  is  the  deformation  part  of  stiffness,  Kyo  is  the  volumetric  part  of  stiff¬ 
ness,  and  Kc  is  the  total  elemental  stiffness  (volumetric  and  deformation  parts 
combined).  For  a  four-noded  quadrilateral  element,  the  volumetric  strain  is  lin¬ 
early  distributed  within  the  element.  The  penalty  function  constraint,  fv  kydV , 
requires  the  value  of  ey  to  be  zero  at  every  point  of  the  element  because  of  the 
‘squared’  term.  This  over-constraint  condition  is  relaxed  by  using  the  one  point 
numerical  integration  scheme,  which  is  applied  at  the  centroid  of  the  element. 

For  identifying  the  tension  and  compression  zones,  the  analysis  is  based  on 
the  volumetric  component  of  the  stiffness  matrix.  It  has  been  observed  that  a 
distinct  separation  exists  between  the  eigenvalues  of  the  shear  deformation  part 
and  the  volumetric  deformation  part  in  the  elemental  level  eigenvalue  represen¬ 
tation,  as  shown  in  Fig.  51.  This  figure  shows  the  real  eigenvalue  spectrum  of 
an  element  during  the  axisymmetric  compression  of  a  disk.  In  general,  the  de¬ 
formation  of  the  material  can  be  characterized  by  three  modes:  the  rigid  body 
mode,  the  shear  mode,  and  the  volumetric  mode.  The  rigid  body  mode  is  due  to 
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Fig  51.  Eigenvalue  spectrum  of  a  generic  elemental  stiffness  stiffness  matrix 


the  rigid  movement  caused  by  adjacent  elements.  The  shear  mode  is  due  to  pure 
deformation  of  the  material.  The  volumetric  mode  is  a  result  of  bulk  deforma¬ 
tion,  which  is  not  practical  to  have  in  any  metal  forming  process.  Therefore,  in 
closed  die  forging  the  volume  is  maintained  constant  during  deformation.  This 
is  achieved  by  using  a  large  penalty  constant  as  explained  earlier.  In  the  elemen¬ 
tal  level  eigenvalue  distribution,  at  least  one  eigenvalue  and  its  corresponding 
eigenmode  is  due  to  the  volumetric  component  of  stiffness.  Because  the  penal¬ 
ty  constant  is  very  large  and  applied  to  the  volumetric  deformation  component, 
the  corresponding  eigenvalue  will  be  the  largest  in  the  real  eigenvalue  spectrum. 
The  eigenvalues  are  calculated  by  first  computing  the  stiffness  matrix  Kyo  (Eq. 
(9.2))  and  using  an  IMSL  subroutine  DEVCSF  to  extract  these  values  from  this 
matrix. 


Once  the  eigenvalue  representing  the  volumetric  deformation  component  is 
identified,  the  corresponding  eigenvector  can  easily  be  calculated.  The  elemental 
level  representation  of  the  eigenvalue  problem  is: 


Kyo<£  = 


\L 


(. KC)dA 


<t> 


(9.3) 


where  K  is  the  penalty  constant,  <f>  is  the  modal  matrix  having  dimensions  corre¬ 
sponding  to  the  number  of  degrees  of  freedom  (8  X  8)  of  each  element,  and  C  is 
the  volumetric  gradient  vector  of  size  1X8.  Eq.  (9.3)  is  of  the  form  Ax  =  Ax, 
which  is  the  general  form  of  a  standard  eigenvalue  problem  where,  A  is  the  eigen¬ 
value  matrix,  and  x  is  the  modal  matrix.  Further,  C  can  be  determined  from 
the  expression  for  the  volumetric  strain  rate  'ey,  as  described  below  [27]: 


—  £z  +  £y  + 


(9.4) 


♦ 


4. 
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In  terms  of  nodal  velocities,  the  volumetric  strain  rate  is  expressed  as, 


•eV  =  CTV  =  Cm  i  =  1,2,3 


(9.5) 


where  C{  =  Bu+ B2 i+ B3i  and  lB ’  refers  to  the  elements  of  the  strain  rate  gradient 
matrix  corresponding  to  the  x,  y  and  z  coordinates.  Now,  the  eigenmode  problem 
may  be  formulated  as, 


Se  =  C<f>\rot 


(9.6) 


where  Se  is  a  scalar  used  for  identifying  the  tension  and  compression  zones,  and 
(j>V  is  the  volumetric  deformation  mode  of  the  elemental  stiffness,  a,  the  scalar 
term,  is  calculated  from  the  volumetric  eigenmodes  as, 


(9.7) 


where  V  is  the  elemental  nodal  velocity  vector.  The  value  of  Se  determines 
whether  the  element  is  in  tension  or  compression.  The  general  rule  is, 

I f  C^>y<£yTV  >  0  the  element  is  in  tension 

If  C </>y<£yTV  <  0  the  element  is  in  compression 

The  product  of  C  and  V  gives  the  volumetric  strain  rate  of  the  element  which  is 
a  scalar,  and  the  product  of  the  eigenvectors  (<£y)  is  another  positive  scalar,  rep¬ 
resenting  the  extent  of  the  deformation.  Because  the  product  of  these  two  scalar 
quantities  gives  the  volumetric  strain  of  the  element,  Se  can  indicate  whether  the 


element  is  in  tension  or  compression. 


Global  level  eigenmode  analysis 


As  explained  in  chapter  2,  the  equilibrium  equations  dealt  with  in  metal 
forming  processes  are  nonlinear,  and  the  solution  to  these  equations  is  obtained 
in  an  iterative  manner.  The  velocity  solution  is  actually  obtained  by  lineariz¬ 
ing  the  equilibrium  equations  and  applying  suitable  convergence  criteria.  The 
linearization  is  achieved  by  a  Taylor  series  expansion  (Eq.  (2.14))  about  the 
assumed  velocity  solution  as  shown  below: 


dir 


L5v/Jv= 


+ 


v=v0 


527T 


A  vj  =  0 


(9.8) 


v=v0 


[ dvjdvj 

where  Vo  is  the  converged  solution  point,  and  A  vj  is  the  first  order  correction  of 
the  velocity  v.  Eq.  (9.8)  can  be  rewritten  by  taking  the  first  derivative  term  to 
the  right  hand  side  as, 

KAv  =  f  (9-9) 


The  linearized  form  of  the  system  of  equations  is  given  as, 

K5V  +  K^Av  =  F  +  f  (9.10) 

where  Ks  is  the  secant  stiffness  matrix  at  the  converged  solution  point  of  the 
previous  iteration,  Ky  is  the  gradient  stiffness  matrix  at  the  current  iterative 
step,  F  is  the  force  vector,  and  f  is  the  change  in  force  due  to  a  change  in  the 
nodal  velocity. 


4 


Using  the  discrete  representation  of  the  quantities  involved  in  the  stiffness  ^ 

equation  developed  in  chapter  2  (Eq.  (2.12)),  we  have  [27], 

d-K  _  dirp  chjp  diTSf  (9  -Q)  * 

dvj  dvj  dvj  dvj 
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where  the  first  term  on  the  right  hand  side  of  the  above  equation  is  the  load 
vector  due  to  deviatoric  stresses,  the  second  term  is  the  load  due  to  hydrostatic 
stress,  and  the  third  term  is  the  applied  nodal  point  force.  These  three  terms 
are  added  up  to  give  the  additional  force  term  f .  The  second  derivative  terms  of 
7 r  are  expressed  as  (Eq.  2.21), 

-/  —  =  f  %-PijdV  +  /  (x-Tjr  -  rPlKVKVMPMjdV  +  f  KCjCidV  (9.12) 
du/uj  JV  e  Jv\ede  eJe  Jv 

The  third  term  in  the  above  equation  is  the  gradient  matrix  due  to  volumetric 
deformation.  This  part  of  the  gradient  stiffness  is  easily  separable  from  the 
other  parts  due  to  its  association  with  the  penalty  parameter,  and  is  used  for 
global  level  eigenvalue  analysis.  The  stiffness  terms  are  first  computed  at  the 
elemental  level  and  then  assembled  into  the  global  stiffness  matrix  (Eq.  2.13). 
This  assembled  stiffness  matrix  has  both  the  shear  and  volumetric  components 
together.  Consequently,  the  real  eigenvalue  spectrum  of  the  stiffness  has  two 
parts.  These  parts  again  are  distinct  and  can  be  easily  identified  and  separated 
because  of  the  penalty  parameter  associated  with  the  volumetric  part. 

In  the  global  level  eigenvalue  analysis,  the  first  step  is  to  identify  the  eigen¬ 
values  related  to  Kvo  &nd  calculate  the  corresponding  eigenvectors  matrix 
The  next  step  in  this  procedure  is  to  compute  a  for  each  volumetric  deforma¬ 
tion  mode,  by  taking  the  product  of  the  eigenvector  and  the  nodal  velocities  as 
explained  in  the  previous  section.  This  computation  results  in  a  vector  given  by, 

a  =  $vTV  (9-13) 
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The  volumetric  gradient  matrix,  C ,  for  each  element  is  assembled  as  a  global 
matrix  CG  of  dimension  equal  to  total  number  of  degrees  of  freedom  of  the 
system.  This  matrix  is  sparsely  populated  and  has  values  only  at  eight  places, 
the  rest  of  the  terms  in  a  given  row  being  zeros.  These  eight  values  correspond 
to  the  degrees  of  freedom  of  the  element  under  consideration.  Now,  the  scalar  S 
is  calculated  for  each  element  as, 

S  =  CG§V«  (9-14) 

Because  CG  is  assembled  row  wise,  only  the  degrees  of  freedom  corresponding 
to  the  element  under  consideration  are  made  active.  Eq.  (9.14)  actually  gives 
the  element  volumetric  strain.  Therefore,  the  sign  of  the  scalar,  S,  provides 
information  whether  the  element  is  in  tension  or  compression. 

Numerical  examples 

The  methodology  developed  is  demonstrated  using  two  numerical  examples 
presented  in  this  section.  The  results  are  validated  by  two  methods:  The  actual 
computation  of  the  volumetric  strain  and  the  g  ratio  method  presented  by  Chen 
et.  al.  [83].  The  first  method  is  simple  and  easy  to  implement  because  the 
volumetric  strain  rate  can  directly  be  calculated  using  Eq.  (9-5).  To  obtain 
the  volumetric  strain,  the  calculated  volumetric  strain  rate  is  multiplied  by  the 
simulation  time  step  (size)  which  is  a  positive  quantity.  In  general,  the  volumetric 
strain  determines  whether  volume  is  gained  or  lost,  and  this  in  turn  can  help 
predict  whether  the  element  is  in  tension  or  compression.  If  the  calculated  value 
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(of  volumetric  strain)  is  greater  than  zero,  the  element  is  in  tension.  Otherwise, 
the  concerned  element  is  in  compression. 


V 


>■ 


> 


The  other  validation  method  used  in  this  work  is  the  g  ratio  method  which 
is  briefly  explained  in  this  section.  The  g  ratio  is  defined  as  [83], 

g  =  ^  (9-15) 

O’ 

where  crm  is  the  mean  or  hydrostatic  stress,  and  ct  is  the  effective  stress.  The 
hydrostatic  stress  in  turn  is  calculated  as, 

0m  =  ^  (*1  +  ^2  +  0-3)  (9-16) 

where  <ti,ct2  and  <73  are  the  principle  stresses. 


The  Dynamic  Material  Model  (DMM)  developed  by  Prasad  et.  al.  [84], 
and  Gegel  [85]  provides  useful  macroscopic  information  for  identifying  favorable 
processing  zones  in  the  workpiece  material  during  deformation.  This  procedure 
models  the  relationship  between  the  constitutive  behavior,  hot  workability,  and 
microstructure  development.  The  modeling  is  based  on  four  criteria:  range  of 
strain  rate  sensitivity,  stability  criteria  relating  the  variations  of  strain  rate  sensi¬ 
tivity  with  loge,  the  lower  limit  on  the  temperature  sensitivity  on  flow  stress,  and 
the  stability  criteria  relating  the  variations  of  the  entropy  with  loge.  According 
to  the  authors,  the  rate  sensitivity  should  lie  between  0  and  1,  and  the  lower 
limit,  on  entropy  should  be  1,  for  stable  deforming  operations.  Further,  for  stable 
processing  conditions,  the  following  relations  should  hold  good: 


dm 

dloge 


<  0 


and 


ds 

dloge 


<  0 


(9.17) 
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where  s  is  the  entropy  of  the  system  and  m  is  the  strain  rate  sensitivity  index. 
By  analogy,  Chen  et.  al.  [83],  considered  the  stress  ratio  parameter  ( g )  as  a 
measure  of  stability,  and  concluded  that, 


If  g  <  0, 


ds 

dloge 


the  element  is  stable 


(9.18) 


and 

ds 

If  R  >  0,  — — r  >  0  the  element  is  not  stable.  (9.19) 

dloge 

Hence,  the  g  ratio  can  be  used  as  a  reference  for  dynamically  determining  favor¬ 
able  processing  conditions  during  deformation.  If  the  g  ratio  value  is  less  than 
i,  the  state  of  stress  is  considered  compressive,  otherwise,  the  state  of  the  stress 
is  considered  tensile. 


Example  1:  H-block  compression 

In  this  example,  the  closed  die  forging  of  a  H-block  is  simulated  using  the 
nonlinear  finite  element  code  ALPID  (Fig.  52).  Isothermal  simulation  conditions 
were  assumed,  with  a  initial  billet  temperature  of  1700°F.  A  constant  shear  fric¬ 
tion  factor  (m=0.3)  was  assumed  at  the  die/billet  interface.  A  Titanium  alloy  was 
used  as  the  workpiece  material.  A  ram  velocity  of  0.12  mm/ sec  was  used  at  the 
start  of  the  simulation.  Because  of  the  symmetry  of  deformation,  only  a  quarter 
section  of  the  workpiece  is  modeled  and  analyzed,  while  maintaining  symmetric 
boundary  conditions  along  the  X  and  Y  axes.  The  simulation  is  carried  out  for  15 
steps,  and  the  billet  deformation  is  predicted.  The  results  of  the  elemental  and 
global  level  analyses  are  compared  with  the  stress  ratio  method,  and  presented 


Fig  53.  Stability  predictions  for  H-Block  compression  (Step  1) 

(a)  Stress  ratio  method  (b)  Elemental  analysis  (c)  Global  analysis 
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Fig  54.  Stability  predictions  for  H-Block  compression  (Step  15) 

)  Stress  ratio  method  (b)  Elemental  analysis  (c)  Global  analysi: 


in  Figs.  53  and  54.  Figs.  53a  and  54a  show  the  predictions  of  the  stress  ratio 
method.  Figs.  53b  and  54b  show  the  results  of  elemental  level  analysis,  while 
Figs.  53c  and  54c  show  the  results  of  global  level  analysis.  The  elemental  level 
results  exactly  match  the  prediction  made  by  the  stress  ratio  method.  In  case 
of  the  global  level  analysis,  the  results  match  in  most  of  the  elements  except  for 
three  or  four  elements  which  are  in  a  state  of  transition  between  the  compression 
and  tension  zones  as  indicated  in  Figs.  53c  and  54c.  This  discrepancy  could  be 
due  to  the  coupling  of  modes  of  deformation  in  the  billet  (i.e.,  coupling  of  the 
deformation  and  shear  modes). 

Example  2:  Disk  forging 

In  this  example,  the  closed  die  forging  of  a  circular  billet  of  radius  52  mm 
is  simulated  under  plane  strain  conditions.  Two  dies  are  used  in  simulating  the 
process.  The  top  die  is  given  an  initial  velocity  of  1.0 mm/ sec,  while  the  bottom 
die  is  kept  stationary.  A  half-model  is  used  for  simulation  of  the  process  because 
of  the  symmetry  in  deformation  (Fig.  55).  Symmetrical  boundary  conditions  are 
applied  for  the  nodes  along  the  Y-axis.  The  billet  is  modeled  with  97  elements. 
A  rate  sensitive  billet  material  having  the  following  constitutive  relationship  is 
used  for  the  simulation: 

a  =  Kentm  +  c  (9.19) 

For  the  current  example  the  constant  K  is  10.0,  the  strain  rate  hardening  sensi¬ 
tivity  index,  m,  is  0.1,  the  strain  hardening  index,  n,  is  0  and  the  constant  c  is 
0.  The  simulation  is  carried  out  under  isothermal  conditions  for  10  steps,  each 


* 


> 


X-fiXIS 


> 


Fig  55.  Plane  strain  compression:  Finite  element  model 
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Fig  56.  Stability  predictions  for  plane  strain  compression  (Step  1)  ^ 

(a)  Stress  ratio  method  (b)  Elemental  analysis  (c)  Global  analysis 
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Fig  57.  Stability  predictions  for  plane  strain  compression  (Step  10) 
(a)  Stress  ratio  method  (b)  Elemental  analysis  (c)  Global  analysis 
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of  step  size  0.1  sec.  The  elemental  level  analysis  results  again  compared  very 
well  with  that  of  the  stress  ratio  method.  The  global  level  eigenmode  analysis 
results  along  with  the  stress  ratio  results  are  presented  in  Figs.  56  and  57.  The 
results  of  the  global  analysis  matched  that  of  the  stress  ratio  method  for  most 
of  the  elements  except  the  ones  indicated  Figs.  56c  and  57c.  The  uncharacter¬ 
istic  behavior  of  some  of  these  elements  is  again  attributed  to  the  coupling  of 
deformation  modes  within  the  workpiece  material. 

Summary 

In  this  work  an  attempt  was  made  to  relate  the  deformation  modes  of  the 
forging  system  with  the  tensile  and  compressive  behavior  of  the  material.  Iden¬ 
tifying  the  tension  and  compression  zones  can  be  directly  related  to  the  stability 
of  the  system  under  deformation.  Because  identifying  the  zones  is  carried  out 
using  the  stiffness  matrices,  the  predictions  depend  closely  on  the  stability  of  the 
stiffness  matrices.  Hence,  the  volumetric  modes  calculated  in  a  particular  step 
can  be  used  for  several  steps  by  calculating  the  stability  margins  of  the  stiffness 
matrices  in  that  step  and  by  confirming  that  in  the  succeeding  steps  the  stiffness 
matrices  do  not  exceed  the  stability  margins. 

The  elemental  analysis  showed  a  clear  separation  of  modes,  namely,  the  rigid 
body  mode,  shear  deformation  mode  and  the  volumetric  deformation  mode.  Each 
one  is  characterized  by  a  different  location  in  the  real  eigenvalue  spectrum.  The 
predictions  made  using  the  elemental  level  analysis  matched  very  well  with  the 


predictions  of  the  stress  ratio  method  and  with  the  volumetric  strain  calculation. 
The  global  level  formulation  was  also  tested  and  compared  with  the  stress  ratio 
method.  There  was  some  discrepancy  in  the  predictions  of  the  two  methods. 
These  differences  may  be  attributed  to  the  coupling  of  shear  and  volumetric 

deformations. 

Out  of  the  three  modes  of  deformation,  the  rigid  body  modes  clearly  separate 
from  the  other  two  modes.  Between  the  shear  and  volumetric  modes,  a  coupling 
exists  for  certain  problems  whose  geometry  and  processes  are  complex.  From 
the  stiffness  calculations,  it  is  clear  that  the  bulk  and  volumetric  strain  rates  are 
involved  in  the  deformation  process.  Both  of  these  are  functions  of  the  strain  rate 
components  ex,  ey,  and  ez.  Due  to  this  nature  of  the  stiffness  matrix  coupling 
exists  between  the  shear  and  the  volumetric  deformation  modes.  This  coupling 
is  not  shown  very  well  in  the  elemental  analysis  because  each  element  is  treated 
individually.  Because  the  global  analysis  is  based  on  the  volumetric  deformation 
of  the  entire  bidet,  the  coupling  of  the  modes  is  more  pronounced. 


CHAPTER  10 


•4 


Summary  and  Conclusions 

The  forging  of  complex  geometries  using  expensive  alloys  under  optimum  pro¬ 
cessing  conditions  is  critical  to  produce  defect-free  and  cost-effective  products. 
This  work  presents  an  innovative  methodology  for  process  modeling  and  control 
using  a  multidisciplinary  approach  based  on  nonlinear  finite  element  methods 
and  optimal  control  theory.  From  the  finite  element  governing  equations  of  the 
metal  forming  system,  a  state  space  model  is  built.  This  model  only  retains  the 
critical  states  of  the  system  (after  condensing  out  from  the  system  the  unwanted 
and  non-critcal  states),  and  has  lower  number  of  degrees  of  freedom  as  compared 
to  the  corresponding  finite  element  model.  This  model  is  further  reduced  us¬ 
ing  sophisticated  model  reduction  methods  if  required.  The  state  equations  are 
solved  using  the  LQR  design  scheme  to  obtain  the  optimal  process  parameters 
(die  velocity  and  initial  die  temperature). 

This  methodology  has  been  implemented  by  means  of  a  computer  program, 
COPP,  developed  during  the  course  of  this  work.  COPP  is  built  using  ALPID 
subroutines,  and  can  be  used  with  isothermal  and  nonisothermal  forgings  to  ob¬ 
tain  optimal  ram  velocity  schedules  and  temperature  trajectories.  Several  design 
problems  have  been  solved  to  substantiate  and  validate  the  methodolgy  devel¬ 
oped.  Both  strain  rate  and  nodal  temperature  have  been  effectively  controlled 
to  satisfy  the  design  requirements,  and  the  performance  of  the  process  has  been 
considerably  improved.  From  this  study,  the  following  conclusions  can  be  drawn: 
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1.  The  condensed  state  space  model  describes  the  forging  process  behavior  very 
well  as  discussed  in  section  6.1  (model  validation)  of  this  report.  This  is  fur¬ 
ther  supported  by  the  numerical  examples  described  in  sections  6.2  through 

V 

6.5. 

2.  Optimal  die  velocity  profiles  may  be  used  to  control  the  strain  rate  in  the 
workpiece,  as  well  as  reduce  the  temperature  gradient  and  temperature  range 
for  any  given  forging  process.  Typically  5  %  to  10  %  reduction  in  temperature 
gradient  was  observed.  At  the  same  time,  the  temperature  range  reduction 
varied  from  15  %  to  20  %  for  the  examples  considered  in  this  work. 

4.  Optimizing  the  initial  die  temperature  directly  influences  the  die-workpiece 
boundary  temperature  and  reduces  the  temperature  range  in  the  billet.  For 
the  numerical  examples  considered  in  this  work  the  temperature  range  was 
typically  reduced  by  about  15  %  due  to  the  optimal  design  and/or  selection 

of  the  initial  die  temperature. 

5.  In  some  cases,  additional  benefits  such  as  a  reduction  in  die  loads  and  process 
time  were  also  observed.  Depending  on  the  processing  conditions  chosen,  the 
process  time  was  reduced  by  upto  25  %  in  some  cases,  while  the  die  load  in 
most  examples  was  reduced  by  15  %  to  20  %. 

Among  the  model  reduction  methods  studied,  the  balanced  model  reduction 
(BMR)  method  was  found  to  be  the  most  effective  in  representing  metal  forming 
processes.  In  future,  other  order  reduction  methods  also  need  to  be  explored  for 
*  this  purpose.  Some  of  the  issues  that  require  further  investigation  are: 

t 
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i.  Issues  pertaining  to  the  stability  and  performance  robustness  of  the  metal 
forming  system. 

ii.  Extending  this  design  procedure  to  suit  on-line  control  applications. 

iii.  Extending  this  methodology  and  making  it  applicable  to  other  unit  forming 
processes. 

iv  Integrating  preform  shape  design  and  staging  criteria  with  the  process  control 
concepts  developed  in  this  work. 


« 
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Appendix  A  -  Project  Overview 


Metal  forming  processes  generally  involve  the  plastic  deformation  of  a  material 
having  a  relatively  simple  geometry  into  a  product  of  relatively  complex  config¬ 
uration  in  one  or  more  operations.  In  metal  forming  technology,  proper  design 
and  control  of  the  process  requires,  among  other  things,  the  determination  of  the 
deformation  mechanics  involved  in  the  process.  For  example,  in  a  forming  process 
such  as  forging,  the  properties,  quality,  and  integrity  of  the  final  product  are  de¬ 
termined  by  factors  such  as  workpiece/die  geometry,  material  properties,  frictional 
conditions,  ram  velocity  schedule,  and  die/billet  temperature  profiles.  Without 
the  knowledge  of  the  influence  of  these  variables  on  the  process  mechanics  it  would 
not  be  possible  to  design  the  dies  and  equipment  adequately,  and/or  predict  and 
prevent  the  occurence  of  defects. 

The  entire  design  process  is  depicted  in  Fig.  Al.  In  this  flowchart,  the  high¬ 
lighted  box  shows  the  control  methodology  developed  and  presented  in  this  report. 
Using  this  methodology,  optimal  ram  velocity  profiles  and  initial  die  temperatures 
may  be  designed  for  any  generic  forming  process,  with  constraints  on  elemental 
strain  rates  and  nodal  temperatures  within  the  deforming  workpiece.  The  exact 
design  procedure  is  decribed  in  detail  elsewhere  in  this  report.  Fig.  A2  briefly 
describes  the  design  procedure  by  means  of  a  flow  diagram  and  also  indicates  the 
chapter  (of  this  report)  in  which  each  of  the  design  procedures  is  elaborated  upon. 


Fig  Al.  Process  control  of  metal  forming  operations  -  Project  overview 
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Fig  A2.  Description  and  organization  of  the  metal  forming  process  control  work 


Appendix  B  -  User’s  Manual  for  COPP 


A  general  purpose  program  COPP  (Control  of  Optimal  Processing  Parameters) 
has  been  developed  for  designing  optimal  die  velocities  and  initial  die  temperature 
*  for  forging  processes.  This  program  has  been  coded  using  FORTRAN  77,  and  uses 

ALPID  (Analysis  of  Large  Plastic  Incremental  Deformation)  subroutines  for  finite 
element  simulation  and  analysis.  A  control  module  has  been  integrated  into  the 
program  for  implementing  the  optimal  design  strategy  described  in  chapter  2.  Be¬ 
cause  the  control  module  uses  IMSL  subroutines  during  the  design  process,  access 
to  IMSL  subroutines  is  necessary  while  linking  this  program.  This  program  runs 
on  VAX  systems  because  ALPID  makes  use  of  some  VAX  system  commands  dur¬ 
ing  the  simulation.  General  exposure  to  ALPID  and  its  postprocessor  (FEMGRA) 
are  required  for  using  this  program. 

COPP  can  be  used  with  both  isothermal  and  nonisothermal  forgings.  For 
isothermal  forging  operations  COPP  requires  an  ALPID  deformation  input  file,  and 
for  nonisothermal  forgings  it  requires  both  ALPID  deformation  and  temperature 
input  files.  For  more  information  regarding  the  ALPID  input  files  the  reader 
is  referred  to  the  ALPID  user’s  manual  [43].  In  addition  to  these  files  COPP 
also  requires  a  small  input  data  file  called  ‘control.dat’,  which  has  information 
regarding  the  design  requirements  and  element(s)/node(s)  to  be  controlled.  After 
►  the  program  is  run,  it  creates  an  output  file  called  ‘output.dat’  which  contains 

the  designed  (optimal)  velocity,  strain  rate  of  the  ‘control’  element,  initial  die 
t  temperature  adjustment,  and  other  relevant  information.  This  manual  presents 
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the  format  for  the  input /output  files  and  also  contains  a  few  sample  data  files  for 
the  reader’s  reference. 


1 

General  format  for  the  input  file  ‘control. dat’ 

4 


ISO 

ICNTL 

NUMD,  IDDIE,  IDPICK,  NUMELE 
VDOT,  SMAX 

NUME,  NUMVI 
NELMZ(NUME) 

EPSLMT(NUME) 

NPVI(NUMVI) 

NUMT 

NTEMP(NUMT) 

GT(NUMT) 

C0EF(4,NUMDIE) 


Note; 


The  variables  with  parentheses  are  arrays,  and  the  quantity  within  the  parentheses 
gives  the  size  of  the  array(s). 


4 
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Description  of  the  input  variables 


ISO 

ICNTL 

NUMD 

IDDIE 

IDPICK 


NUMELE 

VDOT 

SMAX 


O-Isothermal  process,  1— Nonisothermal  process 

0-Simulation  only  (No  design  required) 

1-Optimal  design  required 

(This  means  die  velocity  and/or  initial  die  temperature  need  to  be  designed) 

:  Number  of  control  input(s) 

1- Only  die  velocity  (or  initial  die  temperature)  needs  to  be  designed 

2- Both  die  velocity  and  initial  die  temperature  need  to  be  designed 

:  Number  of  dies  requiring  initial  die  temperature  design 

:  0-User  selects  the  ‘element  of  interest’  and/or  ‘critical  node’ 

1- User  selects  the  ‘critical  node’  and  the  program  selects  the 
‘element  of  interest’ 

(The  program  normally  picks  the  non-perimeter  element  having  the  largest 
strain  rate  in  the  billet  as  the  ‘element  of  interest’) 

2- User  selects  the  ‘element  of  interest’  and  the  program  selects  the  ‘critical  node’ 
(The  program  normally  picks  the  node  having  the  lowest  temperature  in 

the  billet  as  the  ‘critical  node’) 

3-  Program  selects  the  ‘element  of  interest’  and  critical  node 

:  Number  of  elements  used  to  discretize  the  billet 

:  Machine’s  maximum  acceleration  capability 
(This  restricts  the  amount  by  which  the  die  velocity  can  change 
in  a  given  time  step) 

:  Stroke  at  which  the  simulation  is  to  be  stopped 


NUME 


:  Number  of  ‘elements  of  interest’ 


NUMVI 


:  Number  of  nodes  constituting  the  ‘element(s)  of  interest’ 

(If  NUME  is  1,  NUMVI  would  be  4,  for  a  four-noded  quadrilateral  element) 

NELMZ  :  The  element  numbers  of  the  elements  of  interest  separated  by  commas 
(Maximum  of  NUME  sets) 

EPSLMT  :  The  desired  effective  strain  rate  values  for  the  (corresponding)  ‘elements 
of  interest’  separated  by  commas  (Maximum  of  NUME  columns) 


4 


NPVI  :  The  node  numbers  constituting  the  ‘element(s)  of  interest’ 

(separated  by  commas  up  to  NUMVI  columns) 

(Note:  The  node  numbers  for  any  element  have  to  be  in  the  ccw  direction) 
NUMT  :  Number  of  critical  nodal  temperatures  to  be  controlled 


NTEMP  :  The  node  numbers  of  the  ‘critical’  nodes, 

(separated  by  commas  up  to  NUMT  columns) 

GT  :  Desired  temperature  gradient  reduction(s)  (<  0.03) 

(separated  by  commas  up  to  NUMT  columns) 

COEF  :  Coefficients  of  the  desired  temperature  profile 

(4  coefficients  are  needed  for  a  3rd  order  polynomial.  These  have 
to  be  specified  in  sequence  starting  with  the  constant  term 


Note: 


1.  In  the  above  nomenclature  ‘element  of  interest’  refers  to  the  element  whose 
strain  rate  is  being  controlled  (or  needs  to  be  controlled).  Perimeter  elements 


# 

i 
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are  generally  not  picked  as  ‘elements  of  interest’.  Also,  ‘critical  node’  refers 
to  the  die-contacting  boundary  node  whose  temperature  is  being  tracked  (or 
^  needs  to  be  tracked).  The  present  version  of  the  program  can  have,  at  most, 

one  ‘element  of  interest’  and  one  ‘critical  node’  during  the  design  process. 

k 

2.  IfISO=0, 

All  temperature  related  inputs  may  be  omitted. 

3.  If  ICNTL=0, 

EPSLMT,  GT  may  be  omitted. 

4.  If  IDDIE=0, 

COEF  may  be  omitted. 

5.  If  IDPICK=1  or  3, 

NUME,  NUMVI,  NELMZ,  NPVI  may  be  omitted. 

6.  If  IDPICK=2  or  3, 

NUMT,  NTEMP  may  be  omitted. 

* 
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Sample  data  file:  Isothermal  case 
(Die  velocity  design  only) 


ISO 

0 

ICNTL 

1 

NUMD, 

1, 

IDDIE,  IDPICK,  NUMELE 

0,  0,  50 

a 

VDOT, 

15, 

SMAX 

0.5 

NUME, 

1, 

NUMVI 

4 

NELMZ(NUME) 

23 

EPSLMT(NUME) 

0.5 

NPVI(NUMVI) 

27,28,33,34 


In  the  above  input  file,  element  23  has  been  specified  as  the  ‘element  of  interest’, 
and  it’s  strain  rate  must  be  maintained  at  0.5  1/sec.  Refer  section  6.3  for  more 
information  regarding  the  design  of  die  velocities  for  isothermal  forging  processes. 


i 
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Sample  output  file  (‘output. dat^i  Isothermal  case 
(Die  velocity  design  for  5  simulation  steps) 


*********************************** 

*  PROCESS  PARAMETERS 
*********************************** 


ST*  100 

T*10 

VRAM  EPS 

23  EPS 

0  EPS 

0  EPS 

0  EPS 

0.100 

0.0100 

1.0000 

0.49341 

0.00000 

0.00000 

0.00000 

0.00000 

0.200 

0.0200 

1.0000 

0.46508 

0.00000 

0.00000 

0.00000 

0.00000 

0.549 

0.0533 

1.0481 

0.49065 

0.00000 

0.00000 

0.00000 

0.00000 

0.916 

0.0883 

1.0481 

0.48903 

0.00000 

0.00000 

0.00000 

0.00000 

1.300 

0.1232 

1.1005 

0.51541 

0.00000 

0.00000 

0.00000 

0.00000 

In  the  above  output  file,  the  first  column  gives  the  stroke  (ST),  the  second  column 
gives  the  time  (T),  and  the  third  column  gives  the  optimal  die  velocity  (VRAM) 
corresponding  to  that  stroke  or  time.  It  may  be  noted  that  the  stroke  and  time 
have  been  multiplied  by  suitable  scaling  factors  for  convenience.  The  fourth  column 
in  the  output  file  gives  the  strain  rate  of  element  23,  which  was  specified  as  the 
‘control  element’  in  the  input  data  file.  It  is  also  observed  that  the  strain  rate 
of  this  element  is  maintained  more  or  less  at  0.5  1/sec  as  specified  by  the  design 
requirement.  The  last  four  columns  in  the  output  file  are  meant  for  the  strain 
rates  of  other  ‘control  elements’.  But,  at  present  these  columns  do  not  have  any 
significance  because  the  program  (COPP)  is  capable  of  having  only  one  ‘control’ 

element  at  any  given  time. 


Sample  data  file;  Non-isothermal  case 
(Die  velocity  design  only) 


ISO 

1 

ICNTL 

1 

NUMD,  IDDIE,  IDPICK,  NUMELE 

1,  0,  0,  50 

VDOT,  SMAX 

15,  0.5 

NUME,  NUMVI 

1,  4 

NELMZ(NUME) 

23 

EPSLMT(NUME) 

0.7 

NPVI(NUMVI) 

27,28,33,34 

NUMT 

1 

NTEMP(NUMT) 

18 

GT(NUMT) 

0.01 


A 


In  the  above  input  file,  element  23  is  specified  as  the  ‘control’  element,  and  node 
18  is  specified  as  the  ‘critical’  node.  Also,  it  is  required  to  maintain  the  strain  rate 
of  element  23  at  0.7  1/sec,  and  reduce  the  (nodal)  temperature  gradient  of  node 
18  by  1  %  at  every  time  step.  Refer  sections  6.2  and  6.4  for  more  information 
regarding  the  design  of  die  velocities  for  nonisothermal  forging  processes. 
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Sample  output,  file  (‘output. daV):  Non-isothermal  case 
(Die  velocity  for  5  simulation  steps) 


r 

x 


* 


*********************************** 


*  PROCESS  PARAMETERS 


*********************************** 


ST*  100 
18  TEMP  0 

T*10 

TEMP 

VRAM  EPS  23  EPS  0  EPS  0  EPS  0  EPS  O' 
0  TEMP  0  TEMP  0 

TEMP 

0.00 

0.100 

0.00 

0.0100 

0.00 

1.0000  0.50017 

0.00000 

0.00000 

0.00000 

0.00000 

1734.34 

0.00 

0.00 

0.200 

0.00 

0.0200 

0.00 

1.0000  0.48265 

0.00000 

0.00000 

0.00000 

0.00000 

1733.68 

0.00 

0.00 

0.540 

0.00 

0.0533 

0.00 

1.0193  0.49242 

0.00000 

0.00000 

0.00000 

0.00000 

1731.50 

0.00 

0.00 

0.896 

0.00 

0.0873 

0.00 

1.0495  0.50697 

0.00000 

0.00000 

0.00000 

0.00000 

1729.28 

0.00 

0.00 

1.277 

0.00 

0.1223 

0.00 

1.0895  0.52627 

0.00000 

0.00000 

0.00000 

0.00000 

1727.02 

0.00 

The  output  file  shown  above  has  14  columns.  The  first  two  columns  contain  the 
stroke  (ST)  and  time  (T),  respectivly.  The  third  column  gives  the  optimal  ram 
velocity  (VRAM)  for  the  process,  and  the  fourth  column  contains  the  strain  rate 
of  elment  23  (‘control  element’)  due  to  the  designed  velocity  schedule.  Column  9 
give  temperature  of  node  18  (‘critical  node)  at  every  time  step  of  the  simulation. 
The  columns  containing  zeros  are  not  relevant  at  this  stage  and  may  be  ignored. 
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Sample  data  file:  Non-isothermal  case 
(Die  velocity  design  &  initial  die  temperature  design) 


ISO  ~ 

1  1 

ICNTL 

1 

X 

numd,  iddie,  idpick,  numele 

2,  1,  o,  50 

VDOT,  SMAX 

15,  0.5 

NUME,  NUMVI 

1,  4 

NELMZ(NUME) 

23 

EPSLMT(NUME) 

0.7 

NPVI(NUMVI) 

27,28,33,34 

NUMT 

1 

NTEMP(NUMT) 

18 

GT(NUMT) 

0.01 

COEF(4,NUMDIE) 

1735,-550,355,-25 


The  above  input  file  is  similar  to  the  one  shown  earlier  for  nonisothermal  die 
velocity  design.  But  here,  in  addition  to  die  velocity,  it  is  specified  that  initial  die 
temperature  also  needs  to  be  designed.  In  this  case  there  are  two  control  inputs 
(instead  of  one,  as  in  the  previous  case),  namely  the  die  velocity  and  the  initial 
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die  temperature  adjustment.  Here  again,  element  23  and  node  18  are  specified  as 
the  ‘control’  element  and  ‘critical’  node,  respectively.  The  design  requirements  are 
to  maintain  a  strain  rate  of  0.7  1/  sec  in  element  23,  and  reduce  the  temperature 
gradient  at  node  18  by  1%  at  every  simulation  time  step.  The  last  line  in  the  data 
file  gives  the  coefficients  of  a  third  order  polynomial  (in  sequence)  starting  with 
the  coefficient  of  the  constant  term.  This  polynomial  specifies  the  temperature 
trajectory  to  be  tracked  by  node  18  in  order  to  meet  the  design  goals.  The  reader 
is  referred  to  section  6.5  of  this  report  for  more  details  on  the  initial  die  temperature 
design  procedure. 


Sample  output  file  (‘output.dat’l:  Non-isothermal  case 
(Die  velocity  k  initial  die  temperature  design  for  5  simulation  steps) 


*********************************** 


*  PROCESS  PARAMETERS 


*********************************** 


ST*100 

T*10 

DIE-TEMP+ 

0.100 

0.0100 

0.00 

0.00 

0.200 

0.0200 

0.00 

0.00 

0.550 

0.0533 

304.58 

0.00 

0.936 

0.0883 

182.50 

0.00 

1.382 

0.1251 

179.66 

0.00 

In  the  output  file  shown  above,  the  first  column  gives  the  stroke  (ST),  and  the  sec¬ 
ond  column  gives  the  time  (T).  The  third  column  gives  the  initial  die  temperature 
adjustment  value.  The  average  of  the  initial  die  temperature  adjustments  for  all 
the  time  steps,  needs  to  be  added  to  the  original  initial  die  temperature,  to  get 
the  optimal  initial  die  temperature.  In  addition  to  designing  the  initial  die  tem¬ 
perature,  die  velocity  design  is  also  done  in  this  example.  The  reader  is  referred 
to  section  5.2  and  section  6.5  for  further  information  on  initial  die  temperature 
design. 


Appendix  C  -  User’s  Manual  for  Reduced  Order  Models 


The  reduced  order  state  space  models  described  earlier  are  obtained  using 
the  BMR  (Balanced  Model  Reduction)  method.  A  FORTRAN  subroutine  RED- 
MOD  has  been  developed  for  this  purpose.  The  subroutine  REDMOD  takes  in 
the  full  size  state  space  matrices  and  returns  the  reduced  order  matrices.  This 
subroutine  has  been  integrated  into  the  program  COPP  for  designing  optimal 
process  parameters  using  reduced  order  state  space  models.  The  rest  of  the  pro¬ 
gram  including  the  control  design  strategy  is  similar  to  that  used  in  COPP  for 
the  full-size  state  space  models.  The  program  with  the  capability  to  reduce  the 
size  of  the  state  space  system  has  been  named  MODRED  (Model  Reduction). 

The  following  common  blocks  and  subroutine  call  statements  are  needed  to 
link  the  reduced  order  model  routine  (REDMOD)  to  the  full  state  space  model 
program. 


c 


C  COMMON  BLOCK  FOR  FULL  STATE  MODEL 


c 

COMMON  /FULL/  AA(NN,NN),BB(NN),CC(NN,NN),WW(NN),NDD 
COMMON  /OPTL/  VGOOD(5),B(MS,5),NUMDX 
COMMON  /CSYC/  CSYS(MS,MS),EBI(MS,MS) 

COMMON  /MORE/  IMRE 


" \ 


CALL  REDMOD(IPQ) 


SUBROUTINE  REDMOD(IPQ) 


r 
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C  COMMON  BLOCK  FOR  FULL  STATE  MODEL 


C 

COMMON  /FULL/  AA(NN,NN),BB(NN),CC(NN,NN),WW(NN),NDD 
COMMON  /OPTL/  VGOOD(5),B(MS,5),NUMDX 


a.  In  the  common  block  FULL  (added  in  the  control  subroutine), 

AA  :  A  copy  of  the  Plant  matrix. 

BB  :  A  copy  of  the  Input  matrix. 

CC  :  A  copy  of  the  Output  matrix. 

WW  :  A  copy  of  the  Constant  perturbation  term  in  the  system. 

NN  :  Maximum  allowable  size  of  the  full  system. 

NDD  :  Number  of  rows  in  the  C  matrix  (This  depends  on  the  numer  of  outputs). 

b.  In  the  common  block  CSYC  (added  in  the  control  subroutine), 

CSYS  :  The  reduced  model  of  the  output  matrix  is  written  in  this  array. 


EBI  :  Array  to  store  the  inverse  of  the  transformation  matrix  used  in  the  balanced 
representation. 

c.  In  the  common  block  OPTL  (added  to  the  control  subroutine). 

VGOOD  :  Vector  to  store  the  designed  process  parameter(s). 

B  :  Reduced  order  model  of  the  input  matrix  is  stored  in  this  array. 

NUDX  :  Number  of  design  variables. 

d.  In  the  common  block  MORE  (added  to  the  control  subroutine). 

IMRE  :  If  IMRE  =  1,  model  reduction  is  performed.  Otherwise,  model  reduction  is 
not  performed. 

e.  From  the  CALL  statement. 

IPQ  :  Goes  in  to  the  subroutine  MODRED  as  the  size  of  the  full  state  model  and 
returns  as  the  size  of  the  reduced  model. 

NOTE  :  —  After  the  call  statement  for  generating  the  reduced  order  models,  the 
full  state  vector  and  the  input  matrix  need  to  be  multiplied  with  the  inverse  of 
the  transformation  matrix  EBI. 

The  data  files  required  for  running  MODRED  are  the  same  as  those  used 
for  COPP.  The  only  difference  in  this  case  is  that  the  file  ‘control.dat’  has  an 
additional  input  variable  called  IMRE  appended  at  the  end  of  the  data  file.  This 


variable  acts  as  a  switch  to  turn  the  model  reduction  capability  ‘or 
or  ‘off’  (IMRE  =  0).  The  output  of  this  program  (‘output.dat’)  is 

as  in  the  case  of  COPP  (Appendix  B). 


’  (IMRE  =  1) 
also  the  same 


